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ABSTRACT 

A  Point  Defense  Missile  Simulation  has  been  developed.   This  report  describes 
the  concept  of  such  a  missile,  the  basic  features  of  the  simulation  program 
including  the  integration  routine  and  the  jet  reaction  controllers,  and  pro- 
vides a  FORTRAN  coded  source  program. 
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c.        Corrector  for  state  variables  at  time  step  i 

C, ,       Moment  coefficients  of  baseline  missile 
expressed  in  pitch  axes 
in' 

n' 


C  ,       Force  coefficients  of  baseline  missile 
x       expressed  in  pitch  axes 

y' 

z' 


d  Missile  diameter,   ft 

e   -e_  Quaterion  parameters 

h  Step  size  for  fast  integration  speed,   sees 
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K   ,  K  Moment  amplication  factors   due  to  JRC's 
m       n 

K   ,  K  Force  amplication  factors   due  to  JRC's 
y       z 

m  Missile  mass,   slugs 

m.  Modified  value  of  state  variables  at  time  step     i 

m' .  Derivatives   of  state- variables  evaluated  with     m. 
1  l 

p,q,r  Missile  angular  rates  about  body  axes,   rad/sec 

p.  Predicted  value  of  state  variables  at  time  step     i 
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S  Reference  area,   ft 

T11-T33  Elements  of  transformation  matrix  from  inertial  axes  to  body  axes 

T  Thrust  of  rocket,   lbs 
r 

T__  Thrust  of  JRC's,   lbs 
JRC 

u,v,w  Components   of  missile's   inertial  velocity  in  body  axes,   ft/sec 


in 


Um,v_,w_  Components  of  total  relative  wind  velocity  in  body  axes, 
ft/sec 

u    ,v   ,w  Components   of  true  wind  velocity  in  body  axes,   ft/sec 
www  J  J  )         / 

V  Magnitude  of  relative  wind  velocity,   ft/sec 
x,y,z  Body  axes   coordinate   system 

x      ,y      ,z  Target  position  relative  to  missile   in  body  axes,   ft 

x',y',z'  Pitch  axes   coordinate   system 

X,Y,Z  Inert ial  axes   coordinate   system 

X    ,Y    ,HT  North-south,   east-west,   and  height   inertial  position  of 
missile,   ft 

X^,Y„,Z  Total  forces   along  missile  body  axes,   lbs 

X„,,Y,,,Z„  Total  moments  along  missile  body  axes,   lbs-ft 

M  M  M 

XT  ,YT  ,HTT  North-south,  east-west,  and  height  inertial  position  of 

1   l  target,  ft 

x  Position  of  missile  center  of  gravity,  ft 
eg 

x  Position  of  missile  JRC  thrusters ,  ft 
JKC 

y.,  z.  Predicted  values  of  state  variables  at  time  step  i 

a  Missile  angle  of  attack,  degrees 

e  Quaterion  orthogonality  error 

e_  Angular  error  between  gyro  spin  axis  and  target  line  of 
sight,  rad 

e r  Angular  error  between  coil  housing  axis  and  target  line  of 
s  ight ,  rad 

9  ,$  ,y  Standard  missile  Euler  angles  ,  deg 

9,  Angular  error  associated  with  seeker  gimbal  ring,  rad 

aA  ,ov,  Target  elevation  and  azimuth  angles  relative  to  missile  axes 
A  B 

$  Seeker  roll  angle,  rad 

$  Orientation  of  relative  wind  with  respect  to  thruster  axes,  deg 
J 

$  Orientation  of  relative  wind  with  respect  to  missile  axes,  deg 
w 

Y  Seeker  coil  housing  look  angle,  rad 
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y  Angular  error  associated  with  seeker  gimbal  ring,  rad 

Q  Angular  rotation  rate  of  coil  housing,  rad/sec 

Q  Angular  rotation  rate  of  gyro  spin  axis  about  its  y  axis, 

y  rad/sec 

Q  Angular  rotation  rate  of  gyro  spin  axis  about  its  z  axis, 

Z 

Cly  Angular  rotation  rate  of  seeker  yoke,  rad/sec 
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INTRODUCTION 

A  POint  DEfense  Missile  Simulation  (PODEMS)  program  has  been  developed 
and  this  report  offers  a  description  of  its  basic  features,  structure,  and 
requirements.  Although  rather  straight  forward  in  nature,  this  program 
provides  the  basic  framework  from  which  further  simulations  of  increased 
complexity  and  sophistication  can  be  easily  implemented. 

The  concept  of  a  point  defense  missile  as  defined  by  this  effort  can  be 
best  understood  by  analyzing  a  typical  flight.  The  seeker  of  the  missile 
initially  acquires  a  low-altitude,  high  speed,  incoming  target.  The  surface- 
to-air  missile  launches  vertically,  and  then  immediately  performs  a  rapid 
pitch-over  maneuver  toward  the  target  with  a  consequent  altitude  gain  of  less 
than  500  feet.  The  primary  controllers  for  this  phase  of  the  flight  are  two 
pairs  of  diametrically  opposed  jets  (jet  reaction  controllers -JRC )  aligned 
perpendicular  to  both  one  another  and  the  missile  axis  of  symmetry.  Upon 
attainment  of  an  approximately  horizontal  flight  path,  the  primary  controller 
of  the  missile  transfers  to  typical  aerodynamic  surfaces  (CANARDS)  which  then 
guide  the  missile  to  intercept.  The  maneuvers  of  lift-off  and  pitch-over,  for 
which  the  time  frame  is  1.  -  1.5  sees,  are  of  primary  interest  and  therefore 
are  the  object  of  this  simulation. 

The  important  features  of  the  simulation  are:  acceptance  of  a  vertical 
launch  configuration,  implementation  of  JRC  controlled  maneuvers,  a  detailed 
simulation  of  a  large-look  angle  seeker,  and  a  dual  speed  integration  routine. 

The  basis  of  the  simulation  is  a  dual  speed  integration  routine  using 

Hamming's  predictor -modifier-corrector  formulation  for  the  recursion  equations. 

The  user  has  the  option  of  specifying  which  state  variables  are  integrated  with 

the  two  different  steps  size  h  and  h   .  Additionally,  because  of  the  • 

s 

singularities  evident   in  the  Euler  angles,   four  quaterion  parameters   are 


employed  to  uniquely  represent  the  missile  attitude  for  all  possible 
orientations. 

All  of  the  missile  parameters,  including  the  aerodynamic  data  for  both 
the  baseline  missile  and  the  JRC's,  are  listed  within  the  report.  However, 
the  reader  is  cautioned  against  the  assumption  that  a  particular  missile  is 
being  simulated  for  the  data  are  only  representative  of  this  type  of  missile, 

The  intention  of  this  effort  was  to  provide  a  general  basic  structural 
program  capable  of  simulating  a  missile  as  a  rigid  body  with  the  specific 
subroutines  for  the  aerodynamic  data,  rigid  body  parameters,  etc.  to  be 
supplied  by  the  user  as  required. 


COMPUTER  PROGRAM  EXPLANATIONS 

A.  Definitions  of  Coordinate  Systems 

1.  Inert ial  Coordinate  System 

An  inertially  fixed  coordinate  system  (X,Y,Z)  is  attached  to  the 
earth  with  the  origin  at  ground  zero,  the  X  axis  indicating  north,  the 
Y  axis  indicating  east  and  the  Z  axis  indicating  the  local  vertical  (positive 
downward ) . 

2.  Missile  Fixed  Coordinate  System 

A  body  fixed  coordinate  system  (x,y,z)  is  located  with  its  origin  at 
the  missile  center  of  gravity,  the  x  axis  as  the  missile's  axis  of  symmetry 
(positive  pointing  forward),  the  y  axis  rotated  negative  ^5  from  the  right- 
hand  pitch  canard,  and  the  z  axis  rotated  accordingly.  See  Figure  1. 

3.  Pitch  Axis  Coordinate  System 

The  origin  and  the  x'  axis  of  the  pitch  axes  coordinate  system 
(x'jy'jz')  are  coincident  with  their  counterparts  in  the  missile  fixed  axes 
system,  while  the  z'  axis  always  coincides  with  the  projection  of  the 
relative  wind  vector  onto  the  y  ,  z  plane.  The  angle  $   indicates  the 
relative  rotation  of  (x^y^z1)  with  respect  to  (x,y,z).  See  Figure  2. 

B.  Missile  Position  and  Orientation 

The  coordinates  X  ,  Y  ,  HT  locate  the  missile  center  of  gravity  with 
respect  to  the  inertial  coordinate  system  in  the  north-south,  east-west, 
and  height  above  ground  zero  directions  respectively. 

The  orientation  of  the  missile  axes  with  respect  to  the  inertial  system 
is  monitored  using  the  standard  Euler  angles  Y  ,  0  ,  $  (yaw?  pitch,  roll) 
with  the  order  of  rotation  as  given.  The  resulting  transformation  matrix 
from  inertial  coordinates  to  missile  coordinates  is 


X 


z 

L  _J 


Til  T12  T13 
T21  T22  T23 
T31  T32  T33 


where 


Til  =  cos  Y  cos  9 

T12  =  sin  Y  cos  9 

T13  =  -  sin  9 

T21  =  cos  Y  sin  9  sin  $  -  sin  Y  cos  $ 

T22  =  sin  Y  sin  9  sin  §  +  cos  Y  cos  $ 

T23  =  cos  9  sin  $ 

T31  =  cos  Y  sin  9  cos  $  +  sin  Y  sin  $ 

T32  =  sin  Y  sin  9  cos  $  -  cos  Y  sin  $ 

T33  =  cos  9  cos  $ 

In  the  actual  simulation,  the  Euler  angles  are  not  employed  because  of  their 
singularity  at  9  =  +  90.  ;  however  they  are  calculated  and  output ed  to  aid 
the  program  user  in  visualization  of  the  missile  orientation.  The  following 
equalities  define  the  Euler  angles  when  9  ^  +  90.: 

9  =  sin"   (-TL3) 

Y  =  tan"  (T12/T11)     (k   quadrant  tan"  ) 


-1 


-1, 


$  =  tan"   (T23/T33)     (h   quadrant  tan  ) 

However,  when  9  =  +  90.  ,  Y  and  $  are  undefined  and  the  following  con- 
vention is  adopted 

Y  =  0 

$  =  tan"1  (T21,  T3l)     (h   quadrant  tan'  ) 


C.  Quaterions 

To  avoid  the  singularity  of  the  Euler  angles  at  9  =  +  90.  the  quaterion 
system  of  four  coordinates  is  adopted.  The  introduction  of  an  extra  coordinate 
into  the  system  removes  the  singularity  but  requires  the  addition  of  a  con- 
straint equation  on  the  four  parameters . 

2     2 
The  four  coordinates  are  e  ,e1}e„,e^  with  the  constraint  of  e_  +  e.      + 

2     2 
e?  +  ©o  =  !  (orthogonality).  The  elements  of  the  previously  mentioned 

transformation  matrix  are  functions  of  these  coordinates . 


mm     2     2     2     2 
Til  =  eQ  +  ei  -  e2  -  e3 

T12  =  2(eie2  +  eQe3) 

T13  =  2(eie3  -  eQe2) 

T21  =  2(eie2  -  eQe  ) 

2     2     2     2 
T22  =  e^  +  e2  -  e^  -  e^ 

T23  =  2(e2e3  +  e^) 

T31  =  2(6^  +  eQe2) 

T32  =  2(e2e3  -  e^) 

2     2     2     2 
T33  =  eQ  *  e32  -  e/  -  e/ 

The  differential  equations  for  the  quaterion  parameters  as  functions  of 
the  missile  angular  rates  (p,q,r)  are: 


eQ  ■  -  i   (e^  +  e2q  +  e  r) 

\   =  2    (e0p  "  e3(l  +  e2r^ 


e 


e2  =  §  (e3p  +  eQq  -  e^r) 
e3  =  \   (-  e2p  +  e-jq  +  eQr) 

Mechanization  of  the  constraint  equation  is  achieved  "by  defining  an  error 

,     2     2     2     2, 
e  =  1  -  (eQ  +  e1     +  e2  +  e3  ) 

which  is  a  measure  of  the  violation  of  the  constraint  and  applying  a  correction 

factor  to  each  differential  equation  which  reduces  the  error.  With  a  value  of 

-6 
K  =  1  the  equations  remain  correctly  constrained  within   | e |  £  10"  . 

e^  =  -  h   (enp  +  e_q  +  e„r)  +  Ke^e 


=  -  2  (^P  +  e2q  +  e  r)  +  Ke( 


®l  =  2  (eop  "  e3q  +  e2r^  +  Keie 


e2  =  ■§■  (e3p  +  eQq  -  e^)  +  Ke2e 
=  i  (-  e2P  +  ^  +  e0r)   +  Keoe 


e 


The  required  initial  conditions  on  e~  ,  e.  ,  ep  ,  e_  are  given  as 
functions  of  the  initial  Y  ,  9  ,  §  "by 

eQ  =  cos  (y/2)  cos  (e/2)  cos  ($/2)  +  sin  (y/2)  sin  (e/2)  sin  ($/2)  . 

ex  =  cos  (Y/2)  cos  (e/2)  sin  ($/2)  -  sin  (y/2)  sin  (e/2)  cos  ($/2) 

e2  =  cos  (y/2)  sin  (e/2)  cos  ($/2)  +  sin  (y/2)  cos  (e/2)  sin  ($/2) 

e  =  -  cos  (Y/2)  sin  (e/2)  sin  ($/2)  +  sin  (y/2)  cos  (e/2)  cos  ($/2) 


D.  Differential  Equations  for  Rigid  Body 

With  X_  ,  Yp  ,  Z,.,  defined  as  the  total  forces  on  the  missile  expressed 
in  missile  axes  x,y,z  respectively  and  XM  ,  YM  ?  ZM  defined  as  the  total 


M  '  XM 


JM 


moments  about  the  missile  center  of  gravity  expressed  in  the  same  axis  system, 
the  differential  equations  of  motion  are: 

u  =  rv  -  qw  +  x_/m 


v  =  pw   - 

ru  +  Yp/m 

w  =  qu  - 

pv  +  Zp/m 

.       XM 
P=- 

X 

q  =  - 


pr  (I  -  I  )  +  Y„ 

x    y     M 


pq  (I  -  I  )  +  Z.. 
r  =  +      x    y  M 


xn 


Yi 


HT 


T 


-1 


u 


-w 


e0  =  -  \   (e-jP  +  e2q  +  e^r)  +  KeQe 
e2  =  i   (eQp  -  e  q  +  e2r)  +  Ke^ 
e3  =  i   (e3p  +  eQq  -  e^)  +  Ke2e 


1^  =  i   (-  e2p  +  e^  +  eQr)  +  Ke  e 


It  is  assumed  that  the  missile  is  symmetric  about  the  x  axis,  no  cross 

-1    T 
products  of  inertia  exist  and  I  =  I  .  T   =  T   where  T  is  the  trans - 

Y    z 

formation  matrix  from  (X,Y,Z)  to  (x,y,z). 


E .  Definition  of  the  Relative  Wind  Orientation 

Two  angles  a    and  $   define  the  orientation  of  the  relative  wind  vector 

with  respect  to  the  missile  axis  system  as  shown  in  Figure  3,  where  u „  ,  v  , 

2    2    2  — 
w   are  the  x,y,z  components  of  the  resultant  wind  and  V  =  (u^  +v   +w_  )2 

Along  each  axis  the  resultant  wind  component  is  the  difference  of  the  missile 

inert ial  velocity  and  the  true  wind  for  that  axis. 


Urp  =  u  - 


w 


Vm  =  V  -  V 

T       w 


w_  =  w  -  w 
T       w 

Additionally,   §T  defines  the  orientation  of  the  resultant  wind  with 
respect  to  each  individual  JRC  (Jet  Reaction  Controller).  Positive  $  is 

d 

defined  as  shown  in  Figure  U. 
F.  Aerodynamic  Data 

The  missile  aerodynamics  are  divided  into  two  distinct  categories:   (l)  aero- 
dynamic coefficients  for  the  baseline  missile  (no  JRC)  and  (2)  amplification 
factors  which  represent  the  effect  of  the  JRC  thrusters .  Both  sets  of  data  are 
functionally  dependent  on  a    and  either  $   or  $T  but  not  on  Mach  number. 

The  baseline  aerodynamic  coefficients  Cr  ,Ct  ,  C  ,  ,  C. ,  ,C,  ,C, 

are  given  as  shown  in  Figures  5-6  •  For  the  present  simulation  C  i=  C  ,  =  C*t= 

C  ,  =  0  . 
n' 

The  effects  of  the  JRC  thrusters  are  summarized  by  amplification  factors 
as  follows  (  y  force  and  moment  amplification  factors  are  used  as  examples ) 

YFl       "YF| 

'JRC  on     'JRC  off 

iV    —  m 

y  JRC 

ymI      "  ymI 

R  _   MJRC  on      'JRC  off 
TJRC  (XJRC} 


The  specific  values  of  K  -  K  ,  K  ,  K   as  programmed  are  shown  in 

y  '  z  '  m  '  n       ° 

Figures  7-8.  For  this  simulation  K  =  K  =  0  .  Additionally,  the  effects 
of  the  JRC  jets  are  assumed  to  be  independent. 
G.  Integration  Routine-^ 

Hamming's  predictor,  modifier,  corrector  set  of  recursion  equations  are 
used  for  the  dual  speed  numerical  integration  of  the  problem  state  variables 
The  following  is  a  brief  explanation  of  the  equations . 

For  a  system  of  n  ordinary  differential  equations 

y'  =  f(x,y) 
where 

y'  =  dy/dx  , 

a  sequence  of  the  solution  variables 

yi  =  y  U±)         i  =  1,  2,  ... 

can  be  expressed  as  a  function  of  previous  y.  and  y. ,  .  With 

h  =  x.  .,  -  x.  Hammings  method  is: 
l+l    l 

PREDICT:         pi+1  =  y±^  +  y  (2y^  -  y^  +  2y|_2) 

112 
MODIFY:  »    =  p    -  g£  (p  -  c.) 


»i+1  =  t   (*i+1  .  »1+1) 


COERECT: 


e1+1  =  |  [*±  -  y±.2  +  3h  (m^ 


+  zy[  -  yi_x)] 


FINAL  VALUE:        jr    =  0    +  jfj  (p    -  e   ) 


Each  advance  of  h  in  the  independent  variable  x  requires  two  evaluations 
of  y'  ,  once  for  the  predictor  and  once  for  the  corrector.  The  method  is 
numerically  stable  with  truncation  errors  to  the  order  of  h   . 

The  values  of  y  and  y'  from  the  past  three  intervals  are  necessary, 
thus  a  starting  technique  is  required.  The  conventional  application  of  a 
Uth  order  Runge-Kutta  integration  method  on  the  first  three  steps  was  dis- 
carded in  favor  of  calculating  the  required  state  variables  by  a  Euler  back- 
step.  Specifically,  for  i  =  0  ,- 


y_3  =  y0  -  3hy« 


y_2  =  y0  -  Shy^ 


y-i  =  yo  "  h  yo 


y-2  :=  y-l  =  y0 


This  method  suffers  from  inaccuracy  when  y'   and  y'   differ  appreciably 
from  y'  .  However,  in  this  simulation,  no  variation  in  the  solution  was 
detected  from  the  application  of  the  less  accurate  Euler  backstep  when  com- 
pared with  a  Runge-Kutta  starter. 

An  additional  complexity  was  introduced  by  the  requirement  of  a  dual 
speed  integration  algorithm  because  of  computational  time  considerations. 
Now  there  are  two  systems  of  differential  equations : 

y'  =  f  (x,y,z) 
z»  =  f  (x,y,z) 

with  the  z  equations  requiring  smaller  time  steps  than  the  y  equations 
for  the  same  accuracy  criteria.  With  h  and  h   defined  as  the  smaller 


10 


and  larger  step  sizes  respectively,  figure  9  depicts  the  sequencing  of  the 

algorithm  for^  one  step  of  h   .  A  ratio  of  h  /h  =  5  is  chosen  for  illustra- 

s  s 

tion  although  this  is  variable  at  the  operator's  option. 
H.  Other  Subroutines 
CONTROL  SYSTEM 

The  JRC's  were  assumed  to  he  the  primary  controlling  elements  for  the 
initial  missile  trajectory  and,  therefore,  the  canard  deflection  are  identically 
zero  for  this  phase  of  the  flight. 

Two  control  equations  govern  the  action  of  the  JRC*s,  one  for  each  pair 
of  opposing  jets.  Figure  10  defines  the  jet  numbers  and  orientations.  For 
illustration  the  control  of  jets  1  and  3  is  presented.  A  variable  FRMTZ  is 
defined  as  a  function  of  missile -target  relative  position  and  rates.  The 

exact  specification  for  this  equation  is  the  operators  responsibility.  When 
FRMTZ  >  0  jet  3  is  on  while  jet  1  is  off.   If  FRMTZ  <  0  the  reverse 
is  true,  and  when  FRMTZ  =  0  both  jets  are  off. 

As  an  example  of  a  possible  control  equation  consider 

FRMTZ  =  crA  +  K  aA 

where  a.   is  defined  in  Figure  11. 
SEEKER 

This  program  incorporates  a  simulation  of  a  large  look  angle  version  of  a 
present  day  seeker.  The  simulation  was  supplied  by  the  manufacturer  and  was 
only  slightly  modified  to  interface  correctly.  The  system  description  will 
not  be  discussed  here,  only  the  inputs  and  outputs  of  the  subroutine. 

The  following  information  is  required  by  subroutine  SEEKER:  X_  ,  Y-  ,  HT  , 
XT  ,  YT  ,  HTT  ,  [T]  ,  p  ,  q  ,  r  ,  p  ,  q  ,  f  .   The  subroutine  returns:   $  , 

\   »  Y3  '  %   >   eB  ■  €C  '  *1  '  \   '  ny  '  %   '   QG  '  QG   for  outPuttin§  if 

y    z 

desired. 
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RIGID  BODY  PARAMETERS 

All  rigid  body  parameters  (mass  ,1  ,1  ,  C.G.  position)  are  linearly- 
interpolated  between  the  initial  values  at  lift-off  and  the  final  values  when 
the  missile  thrust  motor  is  expended.  The  instantaneous  position  of  the 

center  of  gravity  x    is  defined  relative  to  the  reference  point  for  the 

eg 

aerodynamic  data  as  in  Figure  12.  The  following  table  indicates  the  parameters 
as  used  in  the  program. 


PARAMETER 

LIFT  OFF 

BURN  OUT 

m 

6.7^2 

Ik  710 

slug 

h 

65-1 

k8.1 

slug-ft' 

I 

X 

.U20 

.2U5 

slug-ft' 

X 

-.321 

.U06 

ft 

eg 


s  =  .13635  ft 

d  =  .U167  ft 

XJRC=2-U3Uft 

TJRC  =  k0°'  lbS 

T  =  3000.  lbs 
r 

1+ 
ATMOSPHERE 

Both  the  density  and  acoustical  velocity  of  air  as  functions  of  height  are 
generated  within  ATMOS.  A  linear  interpolation  of  these  parameters  is  based 
on  data  from  an  ICAO  Standard  Atmosphere  Table  at  heights  of  0.  and  1000.  ft, 

Additionally  values  for  the  X  and  Y  components  of  surface  winds  maybe 
entered  as  constant  or  functions  of  altitude  depending  on  the  operator's  pre- 
ference. 
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TARGET 

Subroutine  target  calculates  the  time  history  trajectory  of  the  target 
as  a  function  of  its  initial  inert ial  position,  constant  inert ial  velocity 
components  and  time. 
THRUST 

Missile  thrust  is  assumed  to  be  a  constant  THR  for  a  duration  of  burn 
TBURN,  after  which  THR  =  0   .  . 
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Yaw  Canard 


y 


Pitch  Canard 
End  View  Looking  Forward 
Figure  1.  Definition  of  Missile  Fixed  Coordinate  System 


Projection  of 
Relative  Wind 


End  View  Looking  Forward 
Figure  2.   Definition  of  Pitch  Axes  Coordinate  System 


Ik 


a 


Resultant  Wind 


a  =   cos 


-1 


(Vv) 


^  =  cos"  (wT/(vT  +  wT  )2) 


If  c*  =  0  ,  5  s  0 


V 


Figure  3«  Relative  Wind  Orientation 
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JRC   Thruster  Axis 


Projection  of 
Relative  Wind 


End  View  Looking  Forward 

Figure  h.     Relative  Wind  Orientation 

With  Respect  to  the  JRC  Thruster  Axis 
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C    ,   vs    o/     for      i     =  0    ,    ttA   ,    tt/2 


o 


—T 
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^   -  0    ,    tt/2 


-r 
20 


1— 
30 


- 1- 
1+0 


50 


—I — 

6o 


a     in  degrees 


Figure   5-     C    ,    vs    a     for      $     =  0    ,    ttA    ,    tt/2 
z  W  '  ' 
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m 


Cm,   vs   a    for      ^  =  0  ,   ttA   >   tt/2 


T 
10 


20  30  Ho 

o  in  degrees 


Figure  6.     C    ,    vs   a     for      5=0,   tt/H    ,    tt/2 
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K 


-2    - 

*J  = 

n/2   .  tt 

-_L 

*j  =  ° 

0 

1 

1 

1 

20 


Uo 


6o 


a  in  degrees 


Figure  7.     K     vs   a     for     $T  =  0   ,   tt/2    ,   n 
z  <j 


K 


m 


a  in  degrees 


Figure  8.      K     vs   a     for     §T 

°  m  J 
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=  0    ,    tt/2    ,    TT    . 


h 

s 


xi       xi+l       xi+2      xi+3      Xi+U      xi+5 

1.  At  x  =  x.   evaluate  z!  ,  y J  and  predict  z.  .  and  V'+s  using  h  and 

h   respectively, 
s 

2.  At  x  =  x. ,_   evaluate  z!  ,_  ,  y! ,_   and  correct  z.  ,   and  y.  _  using 

i+I  l+l    l+l  l+l       i+5 

h  and  h   respectively. 

3.  With  y.    now  fixed,  at  x  =  x.    ,  x  _  ,  ...  evaluate  only  z!   , 

l+p  i+±    i+^  i+_L 

z!  „  ,  ...  and  sequentially  predict  z.  „  ,  z.  _  ,  ...  and  correct 

zi+2  '  zi+3  '  * ' *  until  X  =  Xi+5  ' 

k.     Repeat  steps  (l)  -  (3)  for  successive  increments  of  h 

s 


Figure  9«   Dual  Speed  Integration  for  One 
Large  Step  h 
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End  View  Looking  Forward 


Figure  10.  Definition  of  JRC  Orientation. 
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Target 


aA  =  tan"1  (z^/X^) 


°B 


=   tan"   (y  _/x  m) 
wmr  mT 


Figure  11.   Definition  of  Target  Azimuth 
and  Elevation  Angles. 
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APPENDIX  A 
PROGRAM  LISTING 

The  following  is  a  FORTRAN  listing  of  the  simulation  program  designed  for 
compatibility  with  the  United  Computing  Service,  Inc.  time  sharing  system. 
A  typical  input/output  listing  is  also  included. 
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i-AGL       1  SXDF  03/28/73.     15.35.59 

00100      program  sxdfc input* output* tafei,tape2>tafe3) 

00110         COMMON  YC  9),DYC  9 ) ,X ,H,WI NDXB>WI NDYB, VI NDZBjRMACH 

00120+        ,CPHIV,XFA,YFA,ZFA,THP,XCG,RMASS,FCAN,YCAN,XM,YM>ZM,XXI,YYI 

00130+  >VSND>SPHIV,HTE,GUEjALPHA,PHIW*XMAX,TBURNjRHO>VEL.,FSI,THT> 

00140+       FHI,ZC 14),DZC 14) 

00150         COMMON  NDIM,IPRINT>NSTEP,NPRINT,ICOUNT,IRATE,IRAT10>NDIMF 

00160*  COMMUNICATIONS  BETWEEN  SUBROUTINE  CTLSYS#AROJRC 

00170        COMMON/CTLSYJ/  IJETC4) 

00180*  COMMUNICATIONS  BETWEEN  SUBROUTINES  TRNSMT> TARGET, SEEKER 

00190        COMMON/TRANS/T1 1 > Tl 2* Tl 3* T2 \, T22, T23> T31 * T32,T33 

00200*  COMMUNICATIONS  BETWEEN  SUBROUTINES  CTLSYS* SEEKER 

00210         COMMON/SEEKR/  EPB, SPHI 1 ,CFHI 1 , EFBC 

00220*  EGUIVALENCING  COMMON  MISSILE  TERMINOLOGY  TO  THE  STATE  VARIABLES 

00230*  Y(I)  AND  ZCI>  AND  THEIR  RESPECTIVE  DERIVATIVES  DYCI)  AND  DZCI) 

00240        EQUIVALENCE  CU,YC 1 ) ) , C V, YC 2 ) ) , CW,YC 3) > * CXI,YC4) >#<YI *Y< 5> ) » 

00250+        <HTjY<6))#<F*Y(7))j(Q*YC8))*CR*YC9>)# 
00260+       CUD,DYC 1 ) )*<VD*DY(2))*<WD#DYC3> )> 

00270+       CXID,DYC4)).»CYID,DYC5)).»CHTD>DYC6))>CPD,DYC7))>CQD.»DYC8))> 
00280+       CRD,DYC9)) 

00290*  EGUIVALENCING  COMMON  MISSILE  TERMINOLOGY  TO  THE  STATE  VARIABLES 
00300*  Y(I)  AND  Z(I)  AND  THEIR  RESPECTIVE  DERIVATIVES  DY(I)  AND  DZ ( I ) . 
00310  EQUIVALENCE  CWYX,ZC3) )>CPHI 1 ,Zl 4) ) * ( WCZ*Z ( 5 ) ) * CFSI 2>ZC6 ) > » 
00320+  (WGY,Z(7) )>CWGZ*ZC8)),CPSI3,ZC9)),CTHT4,ZC 10) >, 
00330+  CE0,ZC11)),CEWZC12)),CE2,ZC13)  ),CE3>ZC  1  4)  )  »  (PHI  1  D,  DZ  (  4  )  )  , 
00340+  CPSI2D,DZC6)) , CE0D,DZC 1 1 ) ) * C El D> DZ ( 12) ) , C E2D,DZC 13) )* 
0035C+       CE3D,DZC14)) 

00360        DATA  C0N4/57 .2957795 1 /» GRAV/32 .2/ 
00370  DATA  CON1/1.0/ 
00360  PETRI EVE(CTLSYS) 
00390  RETRIEVE(INIT) 
00400  PETRI EVEC THRUST) 
00410  PETRI EVE( AERO) 
00420  PETRI EVEC TRNSMT) 
00430  RETRIEVEC TARGET) 
00440  PETRI EVEC RBFRMT) 
00450  RETRIEVEC INTEG) 
00460  RETRIEVECAROBSE) 
00470  RETRIEVEC SEEKER) 
00480  PETRI EVEC ATMOS) 
00490  PETRI  EVEC  APOJRO 
00500  RETRIEVEC I NTEP2) 
00510  PETRI EVEC OUTFUT1 ) 
00520  REWIND  2 

00530*  INITIALIZE  ALL  PROGRAM  VARIABLES 
00540         CALL  INIT 

00550*  THE  INTEGRATION  ROUTINE  REQUIRES  TWO  EVALUATIONS  OF  THE  STATE 
00560*  VARIABLE  DERIVATIVES  FOR  EACH  TIME  STEF .  ADDI Tl ONALLY, THE  STATE 
00570*  VARIABLES  ARE  SEPARATED  INTO  TWO  CATEGORIES;  ZCI)  FOR  VARIABLES 
00580*  WHICH  UPDATE  AT  EVERY  TIME  STEF  H*  AND  YCI)  FOR  VARIABLES  WHICH 
00590*  UPDATE  ONLY  EVERY  "IRATIO"  TIME  STEPS. 

00600*  CALCULATE  THE  NUMBER  OF  DERIVATIVE  EVALUTI ONS  NECESSARY  FOR 
006  10*  "IRATIO"  NUMBER  OF  TIME  STEPS. 
00620        IRT2=2*IRATI0 
00630      3  DO  2  K=1*IKT2 

00640*  THE  START  OF  THE  DERIVATIVE  EVALUATION  LOOF . 

00650*  CALCULATE  THE  TRANSFORMATION  MATRIX  BETWEEN  THE  INERTIAL  AXIS 
00660*  SYSTEM  AND  THE  MISSILE  BODY  FIXED  AXIS  SYSTEM. 
00670        CALL  TRNSMT 

00680*  CALCULATE  THE  DERIVATIVES  OF  THE  GUATERION  VARIABLES. 
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00690        EFS=1 .+T11+T22+T33-4.*E0*E0 

007  00  E0D=-.5*(E1*P+E2*Q+E3*R>+EPS*E0 
00710  E1D=  .5*(E0*P+E2*R-E3*Q>+EPS*E1 
00720  E2D=  .5*(E0*G+E3*F-E1*R)+EPS*E2 
00730        E3D=  .5*(E0*R+E1*Q-E2*P)+EPS*E3 

00740*  IT  IS  NOT  NECESSARY  TO  REEVALUATE  ALL  OF  THE  DERIVATIVES  WHEN 

00750*  ONLY  THE  Z(I>*  AND  NOT  THE  Y(I>*  VARIABLES  ARE  BEING  UPDATED. 

00760*  "IRATE"  EQUALS  1  IN  THIS  SI TUATI ON*OTHERWI SE  IT  EQUALS  0. 

00770*  WHEN  "IRATE"  EQUALS  1* PORTIONS  OF  THE  LOOP  ARE  BYPASSED  TO 

00780*  PREVENT  RECOMPUTATI ON  OF  NON-UFDATED  VARIABLES. 
00790        IF  (IRATE. EC  .  1  >  GO  TO  20 

00800*  COMPUTE  TARGET  POSITION. 
00810        CALL  TARGET 

00820*  COMPUTE  MISSILE  THRUST. 

008  30        CALL  THRUST 

00840*  COMPUTE  THE  INSTANTANEOUS  VALUES  OF  THE  RIGID  BODY  PARAMETERS. 

008  50         CALL  RBPRMT 

00860*  COMPUTE  THE  PERTINENT  ATMOSPHERIC  PARAMETERS. 

00870        CALL  ATMOS 

00880     20  CONTINUE 

00890*  SIMULATE  THE  SEEKER  DYNAMICS.  THE  SEEKER  VARIABLES  MUST  BE 

009  00*  INTEGRATED  WITH  THE  SMALL  TIME  STEP  TO  AVOID  COMPUTATIONAL 

009  10*  INSTABILITIES*  AND  HENCE  ARE  EVALUATED  DURING  EACH  PASS  THROUGH 

00920*  THE  LOOP. 

009  30        CALL  SEEKER 

00940*  BYPASS  PORTIONS  OF  THE  LOOP  BASED  ON  "IRATE". 

00950        IF  (IRATE. EG. 1)  GO  TO  30 

00960*  DETERMINE  THE  MISSILE  CONTROL  VARIABLES. 

00970        CALL  CTLSYS 

00980*  ADD  THE  EFFECTS  OF  ATMOSPHERIC  WINDS. 

00990        UT=U-WINDXB 

01000        VT=V-WINDYB 

01010        WT=W-WINDZB 

01020*  CALCULATE  MISSILE  VELOCITY,  MACH  NUMBER*  AND  DYNAMIC  PRESSURE. 

0  1030        VEL2=UT*UT+VT*VT+WT*WT 

01040        VEL=SGRT(VEL2) 

0  1050        VEL1=SGRT(VEL2-UT*UT) 

01060        GUE=.5*RH0*VEL2 

01070         RMACH=VEL/VSND 

01080*  DETERMINE  ALPHA  AND  PHI  OF  THE  WIND.  AT  THE  SINGULARITY  OF 

01090*  ALPHA=0*  PHI  IS  DEFINED  AS  =0. 

01100        PRMT=CONl 

01110         IF  (VEL.GT.0.)  PRMT=UT/VEL 

01120        IF  (ABS(PRMT) .GE.l .0)  PRMT=SIGN( CON1 *FRMT) 

01130        ALPHA=C0N4*AC0S(PRMT) 

0  1140        FHIW=0. 

01150         SPHIW=0. 

01160        CPHIW=1. 

01170        IF  (VEL1.EQ.0.)  GO  TO  10 

01180        SFHIW=VT/VEL1 

01190        CFHIW=WT/VEL1 

01200        IF  (ABS(SPHIW) .GE.l .0)  SPHI W=SIGN( CON1 * SPHIW) 

01210         IF  (ABS(CPHIW) .GE.l .0)  CPHI W=SI GN( CON1 * CPHI W) 

01220        PHIW=C0N4*AC0S(CPHIW) 

01230        IF  (SPHIW. LT.0.)  PHIW=-FHIW 

01240*  DETERMINE  THE  AERODYNAMIC  FORCES  AND  MOMENTS. 

01250     10  CALL  AERO 

0  1260*  SUM  THE  AERO*  THRUST*  AND  GRAVITY  FORCES. 

0  1270        XF=XFA+THR+T13*GRAV*RMASS 
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£1260  YF=YFA+T23*GRAV*RMASS 

0  1290  ZF=ZFA+T33*GRAV*RMASS 

0130G*  EVALUATE  THE  MISSILE  STATE  VARIABLE  DERIVATIVES. 

0  1310  UD=R*V-Q*W+XF/RMASS 

01320  VD=P*W-R*U+YF/RMASS 

01330  WD=Q*U-F*V+ZF/RMASS 

01 340  PD=XM/XXI 

01350  FRMT=F*(XXI-YYI ) 

01360  QD=<-R*PRMT+YM)/YYI 

01370  RD=(G*PRMT+ZM)/YYI 

01380  XID=U*T1 1+V*T21+W*T31 

01390  YID=U*T12+V*T22+W*T32 

01400  HTD=-(U*T13+V*T23+W*T33> 

0  1410  30  CONTINUE 

01420*  DETERMINE  WHETHER  OUTPUT  IS  DESIRED.  MIFRINT"=0  INDICATES 

01430*  THAT  THE  INTEGRATION  ROUTINE  HAS  ONLY  FREDICTED  AND  HAS  NOT 

01440*  CORRECTED*  THEREFORE  THE  OUTPUT  IS  MEANINGLESS.  "IFPINT"=1 

0  1450*  INDICATES  THAT  OUTPUTING  IS  POSSIBLE. 

01460  IF  (IPRINT.EG .0)  GO  TO  6 

0  1470*  DETERMINE  IF  THE  DESIRED  OUTPUT  INTERVAL  IS  SATISFIED. 

01480  IF  <MOD(NSTEF,NFFINT> .GT.0)  GO  TO  6 

0  1490*  CONVERT  THE  CUATEPIONS  TO  EULER  ANGLES  FOR  OUTPUTING.  AT  THE 

01500*  SINGULARITY  OF  THE7A=+-90,  PSI  IS  DEFINED  AS  =0  AND  PHI 

01510*  IS  CALCULATED. 

01520  IF  CABSCT13) .GE.l .0)  Tl 3=SI GNC CON1, Tl 3 > 

01530  THT=C0N4*ASIN(-T13) 

01540  IF  (ABS(THT) .EG. 90.)  GO  TO  4 

01550  FSI=C0N4*ATAN2(T12,T1 1 ) 

01560  PHI=C0N4*ATAN2<T23,T33> 

01570  GO  TO  5 

01580  4  PSI=0. 

01590  FHI=C0N4*ATAN2(T21,T31  > 

0  1600  5  CONTINUE 

0  16  10*  KEEP  TRACK  OF  THE  NUMBER  OF  TIMES  OUTPUTING  IS  PERFORMED. 

01620  ICOUN7=ICOUNT+l 

01630*  CONVERT  RADIANS  TO  DEGREES  FOR  OUTPUTING. 

01640  PSI20=PSI2*C0N4 

01650  THT40=THT4*C0N4 

01660  PHI 10=PHI 1*C0N4 

01670  EPB0=EPB*C0N4 

01680  PSI30=PSI3*C0N4 

01690  EPEC0=EPBC*C0N4 

01700*  WRITE  THE  OUTPUT  VARIABLES 

01710  WRITEC2,)   X*U*V*W#XI *YI  *>1T#P*Q*R*P5I a THTjPHI j 

01720+  ALFHA,PHIW,PHI  1 0,PSI 20, PS  I 30, THT40, EPBO, EPBCO, 

0  17  30+  WYX,WCZ,WGY,WGZ,FHI  1  D,F  SI  2D,  (  I  JETC  I  ),I  =  1,4) 

01740*  DETERMINE  WHETHER  OR  NOT  ALL  VARIABLES  ARE  TO  BE  UPDATED 

01750*  DURING  THE  NEXT  FASS  THROUGH  THE  LOOP.  ALL  VARIABLES  ARE 

01760*  EVALUATED  ON  THE  FIPST  AND  SECOND  PASSES  AND  ONLY  THE  ZCI> 

01770*  VARIABLES  ARE  EVALUATED  ON  ALL  SUBSEQUENT  FASSES  UNTIL  AN 

01780*  "IRATIO"  NUMBER  OF  STEPS  HAVE  BEEN  TAKEN. 

01790  6  IRATE=0 

01800  IF  (K.GT.2)  IRATE=1 

0  18  10*  CALL  THE  INTEGRATION  ROUTINE. 

01820  2  CALL  INTEG 

01830*  STOF  THE  FROGRAM  WHEN  THE  REQUIRED  TIME  HAS  ELAPSED. 

01840  IF  (X.LT.(XMAX+H*IRATIO) )  GO  TO  3 

01850*  CALL  THE  PROGRAM  TO  SORT  AND  SEQUENCE  THE  OUTPUT  VARIABLES  FOh 

01860*  FOR  TTY  COMPATIBILITY. 
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01870  REWIND  2 

01880        CALL  OUTPUT 1 

01890        END 
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10  SUBROUTINE    INIT 

20  COMMON   YC    9),DY<    9 ) ,X,H,WINDXB,UINDYB,WINDZB,RMACH 

30+  ,CPHIW,XFA,YFA,ZFA,THR,XCG,RMASS,FCAN,YCAN,XM,YM,ZM,XXI,YYI 

40+  ,VSND,SFHIW, HTE, CUE, ALPHA, PHIW, XMAX, TBURN, RHO, VEL, PSI ,THT, 

50+  FHI,Z( 14),DZC 14) 

60  COMMON    NDIM,IPEINT,NSTEP,NPRINT,ICOUNT, IRATE, IRATICNDIMF 

70*  COMMUNICATIONS    WITH    SUBROUTINE    TARGET. 

80  COMMON/I NITL/XT0,YT0,HTT0,XTD,YTD,HTTD 

90  EQUIVALENCE    ( E0,Z< 1 1 ) ) , ( El ,Z ( 12) ) , ( E2,Z< 1 3) ) * < E3,ZC 14) ) 

100  100    FORMAT<*INFUT   H,XMAX,NFRI NT* IRATI 0* ) 

110  110    FORMATC*INFUT   XT0,YT0,HTT0 ,XTD,YTD,HTTD* ) 

120*  INITIALIZE    INTEGRATION    ROUTINE    CONTROL   VARIABLES. 

130  IRATE=0 

140  ICOUNT=0 

150  IPRINT=1 

160*  SET    DIMENSIONS    OF    YCI)    AND   Z(I)    RESPECTIVELY.    IF    DIMENSIONS 

170*  ARE    INCREASED   ALSO    INCREASE    THE    STORAGE   LOCATION.   DIMENSIONS 

180*  SUBROUTINE    I NTEG . 

190  NDIM=9 

200  NDIMF=14 

210  PRINT          100 

220*  READ   FROM    TTY    THE    STEP    SIZE,    RUN    TIME,     OUTPUT    INTERVAL, STEP    SIZE 

2  30*  RATIO. 

240  READ            ,    H,XMAX,NFRINT,IRATIO 

250*  READ   FROM    TTY    THE    TARGET   POSITION   AND   VELOCITY. 

260  PRINT          110 

270  READ            ,XT0,YT0,HTT0,XTD,YTD,HTTD 

280*  TIME=0. 

290  NSTEF=0 

300  X=0. 

310*  INITIALIZE   ALL    STATE   VARIABLES    TO    0. 

320  DO    1     I=1,NDIM 

330  Y(I)=0. 

340  1    DY(I)=0. 

350  DO    2    I=1,NDIMF 

360  ZCI)=0. 

370  2    DZ(I)=0. 

380  HTE=6. 

390  TBURN=4.5 

400*  INITIAL    EULER   ANGLES    FOR   VERTICAL    ORIENTATION. 

410  FSI=0. 

420  THT=3. 141592653/2. 

430  PHI=0. 

440*  COMPUTE  INITIAL  GUATERION  VALUES  FROM  EULER  ANGLES. 

450  E0=+COS<PSI/2.)*COS(THT/2.)*COS(PHI/2.)+SIN(PSI/2.)*SIN 

460+  (THT/2.)*SINCPHI/2.) 

470  El=+C0S(FSI/2.)*C0S(THT/2.)*SIN(PHI/2.)-SIN(PSI/2.)*SIN 

480+  (THT/2.)*C0S(FHI/2.) 

490  E2=+C0S(PSI/2.)*SIN(THT/2.)*C0S(PHI/2.)+SIN(PSI/2.)*C0S 

500+  (THT/2.)*SIN(PHI/2.) 

510  E3=-C0S(PSI/2.)*SIN(THT/2.)*SIN(PHI/2.)+SIN(FSI/2.)*C0S 

520+  (THT/2.)*C0S(PHI/2.) 

530  RETURN 

540  END 
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10  SUBROUTINE    TRNSMT 

20  COMMON    YC    9)*DYC     9 ) ,X,H,WI NDXB,WI NDYB*WINDZB.,RMACH 

30+  ,CPHIW*XFA,YFA,ZFA,THR,XCG>RMASS,FCAN,YCAN,XM,YM,ZM,XXI,YYI 

40+  #VSND*SFHIW,HTE*GUE*ALPHA*PHIW#XMAX#TBURN*RHO*VEL#PSI*THT* 

50+  FHI,Z< 14),DZ< 14) 

60  COMMON   NDIM*IPRINT.»NSTEF,NPRINT,ICOUNT*IRATE,IRATIO,NDIMF 

70  COMMON/TRANS/T1 WT12,T13>T21,T22,T23,T31,T32,T33 

80  EQUIVALENCE    ( E0,Z< 1 1 ) ) * ( El ,Z< 12) ) , <E2,Z< 1 3) ) » CE3*Z( 1 4) ) 

9  0*    CALCULATE    THE    ELEMENTS    OF    THE    TRANSFORMATION   MATRIX    FROM    THE 

100*    QUATERNION   VARIABLES. 

110  T11=E0*E0+E1*E1-E2*E2-E3*E3 

120  T12=2.*(E1*E2+E0*E3) 

130  T13=2.*CE1*E3-E0*E2) 

140  T21=T12-4.*E0*E3 

150  T22=T11-2.*(E1*E1-E2*E2) 

160  T23=2.*(E2*E3+E0*E1) 

170  T31=T13+4.*E0*E2 

180  T32=T23-4.*E0*E1 

190  T33=T11-2.*(E1*E1-E3*E3) 

200  RETURN 

210  END 
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10  SUBROUTINE    TARGET 

20  COMMON   Y<    9)>DYC    9  >  ,X,H,WI  NDXB,WI  NDYB.»WI  NDZB,RMACH 

30+  *CPHIW*XFA*YFA*ZFA*THP.jXCG*RMASS#FCAN*YCAN*XM*YM*ZMjXXI*YYI 

40+  >VSNDjSPHIW>HTE,QUE>ALPHA,FHIW.,XMAX,TBURN,RHO,VEL,PSI  *THT> 

50+  FHI,Z< 14>,DZ< 14) 

60  COMMON    NDIM,IFRINT.»NSTEP,NPRINT>ICOUNT, IRATE, IRATIO>NDIMF 

70  COMMON/TARG/XTI>YTI>HTTI 

8  0*    COMMUNICATIONS    WITH    SUBROUTINE    I  NIT 

90  COMMON/INI TL/XT0,YT0,HTT0,XTD*YTD,HTTD 

100  COMMON/TRANS/T1 1 > Tl 2, Tl 3* T2 1 , T22, T23, T3 1 » T32* T33 

1 10  COMMON/SEEKRl/SIGA,SIGB>SIGAD,SIGBD 

120  EQUIVALENCE    (U,Y(  1  )  )  ,  <  V, YC  2)  )  ,  (W., Y(  3)  >  ,  (XI  ,Y<  4)  >  ,  <  YI  , Y(  5)  >, 

130+  (HT>Y(6))>(F>Y(7) )>(Q>Y<8) )«CR#Y(9)) 

140  DATA    C0N1/57. 29577951/ 

150*  COMFUTE    THE    INERTIAL    POSITION    OF    THE    TARGET 

160  XTI=XT0+X*XTD 

170  YTI=YT0+X*YTD 

180  HTTI=HTT0+X*HTTD 

190*  COMFUTE    THE    TARGET    FOSITION    RELATIVE    TO    THE   MISSILE    IN    THE 

200*  INERTIAL    AXIS    SYSTEM. 

210  XIMT=XTI-XI 

220  YIMT=YTI-YI 

2  30  HTMT=HTTI-HT 

240*  TRANSFORM    THIS    INTO    MISSILE    BODY    COORDINATES. 

250  XMT=T1 1*XIMT+T12*YIMT-T13*HTMT 

260  YMT=T21*XIMT+T22*YIMT-T23*HTMT 

270  ZMT=T31*XIMT+T32*YIMT-T33*HTMT 

280*  COMPUTE    THE    MOTION    OF    THE    TARGET   AS   VIEWED   FROM    THE    MISSILE 

290*  AXES    USING    THE    VELOCITIES    OF    THE    TARGET    AND    THE   MISSILE   AND 

300*  THE    ROTATIONAL    RATES    OF    THE   MISSILE. 

310  XMTD=T1 1*XTD+T12*YTD-T13*HTTD-U-G*ZMT+R*YMT 

320  YMTD=T21*XTD+T22*YTD-T23*HTTD-V-R*XMT+P*ZMT 

330  ZMTD=T31*XTD+T32*YTD-T33*HTTD-W-P*YMT+Q*XMT 

340*  COMFUTE    THE   AZIMUTH    AND    ELEVATION   ANGLES    OF    THE    TARGET   AS 

350*  SEEN    BY    THE   MISSILE   AND    THEIR    RESPECTIVE   RATES. 

360  SIGA=ATAN2(ZMT*XMT)*C0N1 

370  SIGB=ATAN2(YMT,XMT)*C0N1 

380  SIGAD=CONl*(XMT*ZMTD-XMTD*ZMT)/(XMT*XMT+ZMT*ZMT) 

390  SIGED=CONl*(XMT*YMTD-XMTD*YMT)/(XMT*XMT+YMT*YMT) 

400  RETURN 

4 1 C  END 
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10  SUBROUTINE    THRUST 

20  COMMON    Y(    9),DY<     9 ) ,X*H,WINDXB,WINDYB*WINDZB,RMACH 

30+  ,CFHIU,XFA.,YFA,ZFA,THR,XCG,RMASS,PCAN,YCAN,XM,YM,ZM,XXI.,YYI 

40  +  >VSND>SPHIW>HTE>  QUE*  ALPHA,  PHIW^XMAX.,  TBURN,RHO.»  VEL,PSI  ,  THT, 

50+  PHI,ZC 14>,DZC 14) 

60  COMMON    NDIM*IFRINT*NSTEP*NPRINT>ICOUNT*IRATE*IRATIO*NDIMF 

70  THR=3000. 

80  IF    (X.GT.TBURN)    THR=0 . 

90  RETURN 

100  END 
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10        SUBROUTINE  RBPRMT 

20         COMMON  Y<  9>,DY<  9 ) ,X,H, VI NDXB, WINDYB,WI NDZB,RMACH 

30+       ,CFHIW,XFA,YFA,ZFA,THR.,XCG,RMASS.»PCAN.,YCAN.,XM.,YM,ZM,XXI,YYI 

40+  ,VSND, SPHIW,HTE,QUE, ALPHA* PHI W,XMAX,TBURN,RHO,VEL,F SI >THT, 

50+  PHI,Z< 14)>DZC 14) 

60  COMMON    NDIM,IPRINT,NSTEF,NFRINT,ICOUNT,IRATE>IRATIO.,NDIMF 

70  DATA    RMASS0/6.742/,YYI0/65.1/.,XXI0/.420/,XCG0/-.321/., 

8  0+  RMASSD/-.4506/,YYID/-3.778/,XXID/-.03889/,XCGD/.16152/ 

90  IF    (X.GT.TBURN)    GO    TO    1 

100*  IF  THE  MAIN  ROCKET  IS  STILL  ON,  THE  MISSILE'S  MASS, MOMENTS 

110*  OF  INERTIA  AND  CG  POSITION  ARE  CALCULATED  BASED  ON  THE  INITIAL 

120*  VALUE  PLUS  AN  AVERAGE  RATE  OF  CHANGE  TIMES  THE  ELAPSED  TIME. 

130        RMASS=RMASS0+X*RMASSD 

140        XXI=XXI0+X*XXID 

150        YYI=YYI0+X*YYID 

160        XCG=XCG0+X*XCGD 

170      1  RETURN 

1 8  0        END 
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10  subroutine  atmos 

20  common  y(   9>*dy<  9 ) *x*h*wi ndxb*windyb*windzb*rmach 

30+  >cphiw*xfa*yfa,zfa*thr*xcg*rmass*pcan*ycan*xm*ym,zm*xxi*yyi 

40+  *vsnd*sphiw*hte* que* alpha* phi w*xmax*tburn*rho*vel*p si *tht 

50+  PHI*Z< 14)*DZ<14> 

60  COMMON  NDIM*IFRINT*NSTEP*NPRINT*ICOUNT*IRATE*IRATIO*NDIMF 

70  COMMON/TRANS /Tl 1  * T12*T1 3* T21 *T22*T23*T31 * T32*T33 

80  EQUIVALENCE  (HT*Y<6>> 

90  RHO=.0023769-HT*. 0000000688 

100        VSND=11 16. 89-HT*. 00384 

110*  SET  THE  INERTIAL  COMPONENTS  OF  THE  WIND. 

120        WINDXI=0. 

130        WINDYI=0. 

140*    TRANSFORM    INTO   MISSILE   AXES. 

1 50  WI NDXB=WI NDXI *T1 1 +WI NDYI *T 1 2 

160  WI NDYB=WI NDXI *T2 1 +WI NDY I *T22 

170  WINDZB=WINDXI*T31+WINDYI*T32 

180  RETURN 

190  END 
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10  SUBROUTINE  SEEKER 

20  COMMON   YC    9 ) * DY C    9 ) *X*H*WI NDXB*WINDYB*WI NDZB*RMACH 

30+  *CFHIW*XFA*YFA*ZFA*THR*XCG*RMASS*PCAN*YCAN*XM*YM*ZM*XXI*YYI 

40+  *VSND*SFHIW*HTE*QUE*ALFHA*PHIW*XMAX*TBURN*RHO*VEL*FSI*THT* 

50+  FHI*ZC 14>*DZC 14) 

60  COMMON  NDIM*IFRINT*NSTEP*NPRINT*ICOUNT*IRATE*IRATIO*NDIMF 

7  0  COMMON/TRANS/T1 1  * Tl 2* Tl 3* T21 * T22*T23* T31* T32* T33 

80  COMMON/TARG/XTI*YTI*HTTI 

90  COMMON/SEEKR/EFB*SPHI  1*CPHI  1*EPBC 

100        EQUIVALENCE  ( EL*ZC 1 ) ) * C ER*ZC 2) ) * ( WYX*ZC 3) ) * (PHI  1 *Z ( 4) ) * 

1 10+       CWCZ*ZC5>>*CFSI2*ZC6>>*CWGY*ZC7))*CWGZ*ZC8>)*CPSI3*ZC9>>* 

120+  CTHT4*ZC 10))*CELD*DZC 1 ) ) * C ERD* DZ( 2 ) ) * CVYXD* DZ ( 3) ) * C PHI  ID* DZ 

1 30+  C  4) ) * ( WCZ D*  DZ ( 5 ) ) * C  PSI 2D*  DZ  <  6  > ) * ( WGYD*  DZ ( 7  > ) * ( VGZD*  DZ ( 8 ) ) * 

140+  CFSI3D*DZC9))*CTHT4D*DZC 10>>*CP*YC7>)*CG*YC8)>*CR*YC9)>* 

150+  C  F D*  DY  C  7  )  )  *  C  G D*  DY  (  8  )  )  *  (  RD*  DY  (  9  )  ) 

160  EQUIVALENCE    (XI  * YC 4) ) * (YI * YC 5 ) ) * CHT*YC 6 ) ) 

170  DATA    INIT/0/*TlFSI/.0125/*T2PSI/.0033/* 

180+  RKPSI/-9  60./*EPSIDB/2  5./*RKBPSI/.046/*RN2/-l . /*RPSI / 1 0 .5/* 

190+  RIPSIB/2.3/*RKTPSI/6.54/*T2PHI/.005/*TlFHI/.025/*RKPHI/300./ 

200+  *EPHIDB/25./*RKBFHI/.199/*RPHI/8.8/*RKTPHI/28.1/*RIPHIB/2.8  4 

210+  /*RKTT/350./*RKPTE/.2/*TGYEB/4./*TGZEB/4./*Fl*F2*F3*F4/10.*2 

220+       .* .5* .5/*FSI2L*PSI3L*THT4L/2.35* .21* .2 1 /* CKSTP2* CKSTP3* CKSTF 

230+       4/800. *400.*400./*RIYX*RIYY*RIYZ/. 25* .20* . 10/*HS*RIGZ*RI GY*R 

240+  KHG/3.5* .009* .009* .005/*RI CX*RI CY*RI CZ/.04* .08* .08/*RIDFSI 

2  50+  /.0007/*RKCU2/.8/*TYX2*TYY2*TYZ2/3*0./*TCZ3*TCY3*TCX3/3*0./* 
260+  TCXE*  TCYE*  TCZE/3*0 . /*  TYXU*  TYYU*  TYZU*  TCXU*  TCYU*  TCZU*  TRXU* 
270+                 TRYU*  TRZU*  THXU*  THYU*  THZU*  TGYU*  TGZU*  TGY5*  TGZ5/1 6*0 . / 

280  DATA    C0Nl/3.141592654/*C0N2/6. 283185308/ 

290  XMTI=XTI-XI 

300  YMTI=YTI-YI 

310  ZMTI=-CHTTI-HT> 

320  RMT=SQRTCXMTI *XMTI +YMTI +YMTI +ZMTI *ZMTI > 

330  RNXI=XMTI/RMT 

340  RNYI=YMTI/RMT 

3  50  RNZI=ZMTI/RMT 

360        RNMX=T1 1*RNXI+T12*RNYI+T13*RNZI 

370        RNMY=T21*RNXI+T22*RNYI+T23*RNZI 

380        RNMZ=T31*RNXI+T32*RNYI+T33*RNZI 

390        IF  CINIT.EQ.l)  GO  TO  1 

400        INIT=1 

410        ISVSZ=0 

420        SFSI2=SQRTCRNMY*RNMY+RNMZ*RNMZ) 

430        PSI2=ATAN2CSPSI2*RNMX) 

440         PHI  1  =  0. 

450  IF    CSPSI2.NE.0.)    PHI  1 =ATAN2CRNMZ*RNMY) 

460  1    CONTINUE 

470  SPHI 1=SINCPHI 1) 

480  CPHI l=COSCPHI 1) 

490  SPSI2=SINCPSI2) 

500  CPSI2=C0SCPSI2) 

510  SFSI3=SINCPSI3) 

520  CPSI3=COSCPSI3) 

530  STHT4=SINCTHT4) 

540  CTHT4=C0SCTHT4) 

5  50  SPSIT=SPSI2*CPSI3+CPSI2*SPSI3 

560  CFSIT=CFSI2*CPSI3-SPSI2*SPSI3 

570  EPBPSI=-RNMX*SPSIT+RNMY*CPHI 1 +  CPSI T+RNMZ*SPHI 1*CPSIT 

580  EFBTHT=-CRNMX*CPSIT*STHT4+RNMY*CCPHI 1 *STHT4*SPSI T-SPHI 1*CTHT 

590+  4>+RNMZ*CCPHI 1*CTHT4+SPHI 1 *SFSI T*STHT4> ) 
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600  RNCY=-RNMX*SPSI2+RNMY*CPSI2*CFHI1+RNMZ*CPSI2*SFHI 1 

610  RNCZ=-RNMY*SPHI 1+RNMZ*CPHI 1 

620  EFBC=SGRT(RNCY*RNCY+RNCZ*RNCZ> 

630  EPB=SQRTCEPBPSI*EPBPSI+    EPBTHT*EPBTHT> 

640  RLAMCG=SGRTCTHT4*THT4+PSI3*PSI3> 

6  50*  CALCULATION    OF    GAIN   COMPENSATION    IN   ROLL   AXIS    DRIVE 
660  IF    (ABS(PSI2).LT..579)    GO    TO    10 

670  VCCOMP=SIGNC 1 .0^PSI2) 

680  GO    TO    30 

690  10    IF    (ABSCPSI2) .LT..174)    GO    TO    20 

700  VCC0MP=SIGN<3.3>PSI2) 

710  GO    TO    30 

720  20  VCCOMP=SIGN< 10.5#FSI2> 

730  30  CONTINUE 

7  40  THRESH=.10 

7  50  IF    (ISUSZ.EG.l)    THRESH=.05 

7  60  IF    (RLAMCG.GT. THRESH)    GO    TO    40 
770  IF    CABSCFSI2) .GE. THRESH    )    GO    TO    40 
780  EPSR=0. 

790  ISWSZ=0 

800  GO    TO    50 

8  10  40    EFSR=-THT4*VCC0MP 
8  20  ISWSZ=1 

8  30  50    CONTINUE 

840*  ROLL   AXIS    DRIVE    TORQUE   MOTOR 

8  50  ERD=<EFSR-ER)/T2PHI 

8  60  EPHID=CT1PHI*ERD+ER)*RKPHI 
870  EFHIDL=BOUNDCEPHIDB,EPHID> 
880  EPHNET=EPHIDL-PHI 1D*RKBFHI 
890  RIPHI=EPHNET/RPHI 

9  00  TYXMOT=RKTPH I *BOUND( RI PHI B, RI PHI ) 
9  10*  OUTER    LOOK   AXIS    DRIVE    TORQUE   MOTOR 
920  EPSL=PSI3 

930  ELD=CEPSL-EL)/T2PSI 

940  EFSID=(T1FSI*ELD+EL)*RKPSI 

9  50  EPSIDL=BOUND<EPSIDB,EPSID> 

9  60  EPSNET=EPSIDL-RKBPSI*RN2*PSI2D 

9  70  RIPSI=EPSNET/RP'SI 

980  TCZMOT=RKTPSI*BOUND(RIPSIB>RIPSI ) 

990*  GYRO  TORQUE  CONTROL  EQUATIONS 

1 000  TGYENL=RKTT*EPBPSI -RKPTE*PSI 3 

1010  TGZENL=-RKTT*EPBTHT+RKPTE*THT4 

1020  TGYE=BOUNDCTGYEB,TGYENL> 

1030  TGZE=BOUND<TGZEB,TGZENL) 

1040*    SEEKER  YOKE    DYNAMICS 

1050  WYY=Q*CFHI 1+R*SPHI1 

1060  WYZ=-Q*SPHI1+R*CPHI1 

1070  VYYD=QD*CPHI 1 +RD*SFHI 1 -PHI 1D*( Q*SPHI 1 -R*CPHI 1) 

1080  WYZD=-QD*SPHI 1+RD*CPHI 1-PHI 1D*<Q*CPHI 1+R*SPHI 1 ) 

1090  TYXF1=FRICT(FUPHI1D) 

1100  TYX1=TYXM0T+TYXF1 

1110  TYX=TYX2+TYXU+TYX1 

1 120  TYY=RIYY*WYYD+WYZ*WYX*(RIYX-RIYZ> 

1 130  TYZ=RIYZ*WYZD+WYX*WYY*CRIYY-RIYX> 

1140  TYY1=TYY-TYY2-TYYU 

1150  TYZ1=TYZ-TYZ2-TYZU 

1160  TMX1=-TYX1 

1 170  TMY1=-CTYY1*CPHI 1-TYZ1*SPHI 1) 

1 180  TMZ1=-(TYY1*SPHI 1+TYZ1*CFHI 1 ) 
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1 190  RICC=RICX*CPHI 1*CPHI 1+RICY*SFHI 1*SFHI 1 

1200  RIYXE=RIYX+RICC 

1210  WYXD=<TYX-WYY*VYZ*(RIYZ-RIYY))/RIYXE 

1220  PHI1D=WYX-F 

1230*    SEEKER    COIL   HOUSING    DYNAMICS 

1240  WCX=WYX*CPSI2+WYY*SPSI2 

1250  WCY=-WYX*SPSI2+WYY*CPSI2 

1260  TCZF2=FRICT(F2,FSI2D> 

1270  TCZS2=ST0PSCCKSTP2,PSI2L.,PSI2> 

1280  TCZ2=(TCZM0T-CRN2-1 . > *RI DFSI *WYZD> +RN2+TCZF2+TCZS2-RKCU2* 

129  0+  PSI2 

1300  TCZ=TCZ2+TCZU+TCZ3+TCZE 

1310  WCZD=CTCZ-V.'CX*WCY*(RICY-RICX>)/<RICZ+RN2*RN2*RIDFSI  ) 

1320  PSI2D=WCZ-WYZ 

1330  TCXLI=RICX*<WYYD*SFSI2-FSI2D*(UYX*SPSI2-vnfY*CPSI2>)+WCY*WCZ 

1340+  *CRICZ-RICY) 

1  350  TCYLI=RICY*<VYYD*CFSI2-FSI2D*CVYX*CPSI2+WYY*SFSI2>  >+WCZ*WCX 

1360+  *CRICX-RICZ> 

1370  TCX2=TCXLI-TCX3-TCXE-TCXU 

1380  TCY2=TCYLI-TCY3-TCYE-TCYU 

139  0  TYX2=-(TCX2*CFSI2-TCY2*SPSI2> 

1400  TYY2=-(TCX2*SPSI2+TCY2*CPSI2) 

1410  TYZ2=-TCZ2 

1420*    GYRO    DYNAMICS 

1430  TGY=TGY5+TGYE+TGYU 

1440  TGZ=TGZ5+TGZE+TGZU 

1450  TGYEFF=TGY-HS*V/GZ-CRKHG/RIGZ)*<TGZ+HS*WGY) 

1460  TGZEFF=TGZ+HS*WGY+(RKHG/RIGY)*(TGY-HS*WGZ) 

1470  WGYD=TGYEFF/(P.IGY  +  RKKG*RKHG/RIGZ) 

1480  WGZD=TGZEFF/(RIGZ+RKHG*RKHG/RIGY) 

149  0*    GYRO   HOUSING    DYNAMICS 

150!  V7HY=WGY 

1510  WHZ=WGZ 

1 520  WRX=WCX*CPSI 3+WCY*SPSI 3 

1530  VRY=-WCX*SPSI3+VCY*CPSI3 

1540  WRZ=(WHZ-WRX*STHT4)/CTHT4 

1550  PSI3D=URZ-WCZ 

1560  THT4D=WHY-WRY 

1570  TRZF3=FRICT(F3jPSI3D) 

1580  TRZS3=STOPS(CKSTP3*PSI3LjPSI3) 

1590  TRZ3=TRZS3+TRZF3 

1600  TRZ4=-(TRZ3+TRZU) 

1610  THYF4=FRICT(F4,THT4D) 

1620  THYS4=ST0FSCCKSTP4,THT4L,THT4> 

1630  THX4=-THXU 

1640  THY4=THYF4+THYS4 

1650  THZ4=-( TRZ4+THX4*STHT4) /CTHT4 

1660  THY5=-(THYU+THY4) 

1670  THZ5=-(THZ4+THZU) 

1680  TGY5=-THY5 

1690  TGZ5=-THZ5 

1700  TRX4=-<THX4*CTHT4+THZ4*STHT4) 

1710  TRY4=-THY4 

1720  TRX3=-(TRX4+TRXU) 

1730  TRY3=-(TRY4+TRYU) 

1740  TCZ3=-TRZ3 

1750  TCX3=-(TRX3*CPSI3-TRY3*SPSI3) 

1760  • TCY3=-CTRX3*SPSI3+TRY3*CPSI3> 

1770  TCYE=-TGYE 
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1780 

TCZE=-TGZE 

1790 

RETURN 

1800 

END 
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1810*  FUNCTION    ROUTINE    FOR   A   LIMITER. 

1820  FUNCTION    BOUNDCXL,X> 

18  30  BOUND=X 

1840  IF    CABS(X) .GE.XL)    BOUND=SI GNCXL*X> 

18  50  RETURN 

1860  END 
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1870*  FUNCTION  ROUTINE  FOR  FRICTION. 

1880        FUNCTION  FRICTCFLjX) 

1890        IF  (X)  1*2*3 

1900  1    FRICT=FL 

1910  RETURN 

1920  2    FRICT=0. 

19  30  RETURN 

19  40  3    FRICT=-FL 

19  50  RETURN 

1960  END 
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1970*  FUNCTION  ROUTINE  FOR  STOPS  WITH  COMPLIANCE. 

1980  FUNCTION  STOPSC SLOPE*XT*X) 

1990  STOPS=0. 

2000  DEL=ABS(X)-XT 

2010  IF    (DEL.GE.0.)     STOPS=SLOPE*SI GNC DEL*X) 

2020  RETURN 

2030  END 
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10        SUBROUTINE  SEEKER 

1810*  FUNCTION  ROUTINE  FOR  A  LIMITER. 

1870*  FUNCTION  ROUTINE  FOR  FRICTION. 

1970*  FUNCTION  ROUTINE  FOR  STOPS  WITH  COMPLIANCE. 
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10  SUBROUTINE    CTLSYS 

20  COMMON    Y<    9)*DY<     9 ) *X,H, WINDXB, WINDYB,WINDZB,RMACH 

30+  >CFHIW*XFA,YFA.,ZFA,THR,XCG,RMASS.»PCAN,YCAN,XM,YM.,ZM,XXI,YYI 

42  +  jVSND*SFHIW*HTEjGUE*ALFHAjFHIW*XMAX#TBURN,RHOjVELjFSI>THT* 

50+  FHI,Z( 14)>DZ< 14) 

60  COMMON    NDIM*IPRINT*NSTEP*NFRINT*ICOUNT*IRATE*IRATIOjNDIMF 

70  COMMON/CTLSYJ/    IJETC4) 

8  0  COMMON/SEEKR/EFBjSFHI  UCFHI  UEFBC 

9  0*    COMMUNICATIONS    WITH    SUBROUTINE    TARGET 

100  COMMON/SEEKRl/SIGA,SIGB#SIGAD*SIGBD 

110  EQUIVALENCE    <HT,Y ( 6 ) ) j ( FSI 2,Z C 6) ) » (PSI 2D* DZ ( 6 ) ) 

115  EQUIVALENCE    (PHI  1 D* DZ ( 4) ) 

120  DATA    TB JRC/ 1.2/ 

130  DATA    C0N1/.18/ 

140*  ENABLE    THE    CONTROL    SYSTEM    AT    THE    ENABLE   HEIGHT. 

150  IF    (HT-HTE)     1*2#2 

160  1     IJRC=0 

170  TSJRC=X 

180  GO    TO    3 

190  2    IJRC=1 

200*  DETERMINE   WHETHER    THE    JRC    BURN    TIME   HAS    BEEN    EXCEDED,    AND    IF 

210*  SO*    DISENABLE    THE    JRC'S. 

220  IF    (CX-TSJRC).GT.TBJRC)     I  JRC  =  0 

2  30  3    IF    (I JRC)    4*4,5 

240  4    I JET< 1)=0 

250  IJET<2)=0 

260  IJET(3)=0 

270  IJETC4)=0 

280  GO    TO    6 

29  0*  DETERMINE    THE    JRC    STATES    BASED    ON    THE    CONTROL   LAW. 

300*  PRMTY    CONTROLS    THE    JETS    ON    THE   Y    AXIS   WHILE    PRMTZ    CONTROLS 

310*  THE    JETS    ON    THE    Z    AXIS.    "I JET< 1 ) "= 1    SIGNIFIES    THAT    JET    1    IS    ON. 

320  5    FRMTY=SIGB+CONl*SIGBD 

330  PRMTZ=SIGA+CONl*SIGAD 

340  IF    CFRMTZ)    7>8*9 

350  7    IJET(1)=1 

360  IJET(3)=0 

370  GO    TO    10 

380  8    IJET(1)=0 

390  IJET<3)=0 

400  GO    TO    10 

410  9    IJET(1)=0 

420  IJETC3)=1 

430  10    IF    (FRMTY)     11*12*13 

440  11     IJET(2)=1 

450  IJET(4)=0 

460  GO    TO    6 

470  12    IJET(2)=0 

480  IJET(4)=0 

490  GO    TO    6 

500  13    IJET(2)=0 

510  IJET(4)=1 

520  6    CONTINUE 

530*  THE    CANARD    DEFLECTIONS    ARE    SET    TO    ZERO 

540  FCAN=0. 

550  YCAN=0. 

560  RETURN 

570  END 
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10         SUBROUTINE  AERO 

20        COMMON  Y(  9>,DY(  9 ) ,X,H,WI NDXB,WINDYB,WI NDZB^RMACH 

30+        ,CPHIW,XFA>YFA,ZFA,THR,XCG,RMASS,FCAN.,YCAN.,XM,YM,ZM,XXI,YYI 

40+       ,VSND,SPHIW,HTE,GUE>ALPHA,PHlW*XMAX,TBURN,RHO,VEL.»FSI>THT, 

50+       PHI,ZC 14>,DZC 14> 

60        COMMON  NDIM,IPRINT>NSTEP,NPRINT>ICOUNT,IRATE>IRATIO,NDIMF 

70*  COMMUNICATION  WITH  AROBSE 

80        COMMON/AROB/XFB*YFBjZFB*XMB*YMB#ZMB 

90*  COMMUNICATION  WITH  AROJRC 

100        COMMON/AROJ/YFJjZFJ*YMJ*ZMJ 

110*  DETERMINE  BASELINE  AERODYNAMIC  FORCES  AND  MOMENTS. 

120        CALL  AROBSE 

130*  DETERMINE  FORCES  AND  MOMENTS  DUE  TO  THE  JRC 

140        CALL  AROJRC 

150*  SUM  THE  FORCES. 

160        XFA=XFB 

170  YFA=YFB+YFJ 

180  ZFA=ZFB+ZFJ 

190*  SUM    THE   MOMENTS   AND    TRANSFER    REFERENCE   POINT    TO    THE 

200*  MISSILE   CG. 

210  XM=XMB 

220  YM=YMB+YMJ+ZFA*XCG 

230  ZM=ZMB+ZMJ-YFA*XCG 

240  RETURN 

250  END 
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10  SUBROUTINE   AROBSE 

20  COMMON    Y(     9>*DY(     9 ) *X*H*U I NDXE* WI NDYB* WI NDZB*RMACH 

30+  *CFHIV*XFA*YFA*ZFA*THR*XCG*RMASS*FCAN*YCAN*XM*YM*ZM*XXI*YYI 

40  +  *VSND*SFHIV*HTE*CUE*ALPHA*PHIV*XMAX*TBURN*RHO*VEL*FSI*THT* 

50+  FHI*ZC 14)*DZC 14) 

60  COMMON    NDIM*IFRINT*NSTEP*NPRINT*ICOUNT*IRATE*IRATIO*NDIMF 

70  COMMON/AER01/DFHI*DALFHA 

8  0  COMMON/AROB/XFB*YFB*ZFB*XMB*YMB*ZMB 

90  DIMENSION    CXB< 25* 3) *CYPB< 25* 3) *CZPB< 25* 3) * CLB< 25* 3) * 

100+  CMFB(25*3)*CNPB(25*3) 

110  DATA    CXB    /7  5*0./ 

120  DATA    CYFB/75*0./ 

130  DATA    CZFB/ 

140+  0.*-2.*-4.0*-6.5*-10.0*-13.0*-17.0*-21 .0* -25 .0* -29 . 5* 

150+  -34.0* -38.0* -42.0* -45.0* -48.0* -51 .0*-52.5*-52.5*-51.0* 

160+  -46.0*-45.0*-45.0*-45.0*-45.0*-45.0* 

170+  0.* -2.* -4.  5* -8.0* -12.0* -16.0* -20.  5* -2 5.0* -29. 5* -34. 5* 

180+  -39.0* -43.0* -47 .0* -50 .0* -52 . 5* -52  .5* -48 .0* -45 .0* -49 .0* 

190+  -51 .0* -52.0* -53. 5* -54. 5* -55. 5* -55.5* 

200+  0.* -2.* -4.0* -6. 5* -10.0* -13.0* -17.0* -2 1.0* -25.0* -29 .5* 

2  10+  -34.0* -38.0* -42.0* -45.0* -48.0* -51  .0* -52 . 5* -52 .5* -5 1 .0* 

220+  -46.0* -45.0* -45.0* -45.0* -45.0* -45.0/ 

230  DATA    CLE    /75*0./ 

240  DATA    CMPB/ 

2  50+  0.*-8.*-16.*-26.*-40.*-60.*-80.*-100.*-124.*-142.* 
260+  -164.* -18 6.* -2 10.* -2 30.* -2 56.* -278.* -290.* -290.* -266.* 
27  0+  -206.*-178.*-166.*-160.*-158.*-158.* 

280+  0.*-8.*-20.*-36.*-56.*-76.*-100.*-128.*-158.*-184.* 

290+  -2 14.* -238.* -264.* -282.* -300.* -286.* -250.* -200.* -2 14.* 

300+  -2 18.* -2 18  .  * -2 18.* -2 18.* -220.* -220.* 

310+  0.*-8 .*- 16.* -26.* -40.* -60.* -80.* -100.* -124.* -142.* 

320+  -164.* -186.* -2 10.* -230.* -256.* -278.* -290.* -290.* -266.* 

3  30+  -206.* -178.* -166.* -160.  *-l 58.* -158./ 
340  DATA    CNFB/75*0./ 

350  DATA    CONl/45./*CON2/2.5/*CON3/90./ 

360  DATA    DIA/.4166666667/*S/. 1363538470/ 

370  FHIWT=AM0D(ABSCPHIV)*C0N3) 

380  PHIWTE=FHIWT/C0N1+1 .0 

39  0  IP=PHIWTE 

400  DPHI=FHIVTE-IP 

410  ALFHAE=ALFHA/C0N2+1 .0 

420  IA=ALFHAE 

430  IF    (IA.GE.25)     IA=24 

440  DALFHA=ALPHAE-IA 

450  CALL    INTER2CCXB    (IA*IF)*CXB    ( I A+ 1  *  IF )  *XXB    (IA*IP+1) 

460+  CXB    ( IA+1*IF+1)*CX    ) 

470  CALL    INTER2(CYFB<IA*IF)*CYPB(IA+1*IF)*CYPB<IA*IF+1 ) 

480+  CYFB( IA+1*IP+1)*CYF) 

490  CALL  INTER2(CZPB(IA*IF)*CZPB(IA+1*IP)*CZFB(IA*IP+1) 

500+  CZFB(IA+1*IP+1)*CZP) 

510  CALL  INTER2CCLB  (IA*IP)*CLB  C I A+ 1  *  IP ) * CLB  (IA*IP+1) 

520+  CLB  (IA+1*IF+1)*CL  ) 

530  CALL  INTER2(CMFB(IA*IF)*CMFB(IA+1*IF)*CMFB(IA*IP+1 ) 

540+  CMFB(IA+1*IF+1)*CMF) 

550  CALL    INTER2<CNFB(IA*IF)*CNFB(IA+1*IP)*CNPB<IA*IF+1 ) 

560+  CNPB< IA+1*IP+1 )*CNF) 

57  0  CY=CYF*CPHIW+CZF*SPHIW 

580  CZ=-CYP*SFHIV+CZP*CPHIU     . 

590  CM  =  CMi-*CI-HIW  +  CNf*SPHIW 
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600  CN=-CMF*SFHIV+CNP*CFHIW 

610  XFB=GUE*S*CX 

620  YFB=GUE*S*CY 

630  ZFB=QUE*S*CZ 

640  XMB=GUE*S*DIA*CL 

650  YMB=GUE*S*DIA*CM 

660  ZMB=£UE*S*DIA*CN 

670  RETURN 

680  END 
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10  SUBROUTINE  AROJRC 

20  COMMON    Y<     9)*DYC     9 ) *X*H* VI NDXB* WINDYB* WI NDZB*RMACH 

30+  *CPHIW*XFA*YFA*ZFA*THR*XCG*RMASS*FCAN*YCAN*XM*YM*ZM*XXI*YYI 

40+  *VSND*SPHlV*HTE*GUE*ALF-HA*PHIW*XMAX*TBURN*RHO*VEL*RSI*THT* 

50+  PHI*Z< 14)*DZ( 14) 

60  COMMON    NDIM*IPRINT*NSTEF*NI-RINT*ICOUNT* IRATE* IRATIO*NDIMF 

7  0  COMMON/AERO  1 /DPH 1 * DALPHA 

80  COMMON/CTLSYJ/    IJETC4) 

90  COMMON/AFOJ/YFJ*ZFJ*YMJ*ZMJ 

100  DIMENSION    RKYB<25*3)*RKZBC25*3>*RKMBC25*3)*RKNB<25*3> 

110  DIMENSION    CPE<4)*SPE<4> 

120  DATA    RKYB/7  5*0./ 

130  DATA   RKZB/ 

140+  -1.0*-1.0*-1.0#-1 •0#-l.0j-«99*-.97i-.92j 

150+  -.88*-. 86*-. 81*-. 76*-. 70*-. 62*-.  58*-. 53* 

160+  -.52,-. 53*-. 58, -.66*-. 72*-.  78*-.  80*-.  82*-. 82* 

170+  50*-1.0/ 

180  DATA    RKMB/ 

190+  1.0*1.0*1.0*1.0*1.0*1.03*1.06*1.10* 

200+  1  .16*  1  .22* 1 .29* 1 .38* 1 .47* 1 .54* 1  .60*1  .63* 

210+  1 .65* 1 .60* 1 .48* 1 .30* 1  .18* 1  .08*1  .04*1  .04*1  .04* 

220+  50*1.0/ 

2  30  DATA    RKNB/7  5*0./ 

240  DATA    T JRC/4 50. /*X JET/2. 434/ 

250  DATA    CONl/90./*CON2/2.5/*CON3/180./ 

260  DATA    FHIE1/0./ 

27  0  DATA    CPE/1 .*0.*-l  . *  0 . /* SPE/0 . *  1  .*0.*-l  ./ 

280  ALPHAE=ALPHA/C0N2+1 .0 

29  0  IA=ALPHAE 

300  IF    (IA.GE.25)     IA=24 

310  DALPHA=ALPHAE-IA 

320  RKY=0. 

330  RKZ=0. 

340  RKM=0. 

350  hKN=0. 

360  DO    1    1=1*4 

370  IF    CIJETCI ) .EG.0)    GO    TO    1 

380  PHI J=PHIW-PHIE1-C0N1*(I-1  .) 

390  CP=CPE(I) 

400  SP=SPE(I) 

410        IF  (ABSCPHIJ) .GT.180.)  FHIJ=FHI J-SI GNC 360 . * FHIJ) 

420         PHI JT=AES(PHI J) 

430  PHI JTE=hHI J1/CON1+1 .0 

440         IF=PHIJTE 

450        IF  CIF.GE.3)  IP=2 

460         DPHI=PHI JTE-IP 

470         CALL  INTEH2(RKYB(IA*IF ) *RKYB< I A+l *  IF > *RKYE< I  A* I F  + 1 )* 

480+       RKYB(IA+1*IF+1 )*RKYJ) 

490         CALL  INTER2(RKZB(IA*IP)*RKZB(IA+1*IF)*RKZB(IA*IP+1 )* 

500+       RKZBCIA+1*IP+1 )*RKZJ) 

510         CALL  INTER2(RKMB(IA*IP)*RKMB<IA+1*IP)*RKMB(IA*IF+1 )* 

520+        RKMB(IA+l*Ii-  +  l )*HKMJ) 

530         CALL  INTER2(RKNB<IA*IP)*RKNB(IA+1*IP)*RKNB(IA*IF+1 >* 

540+        RKNBC IA+1*IP+1 ) *RKNJ> 

550        IF  (PHIJ.GE.0.)  GO  TO  2 

560  RKYJ=-RKYJ 

565  RKNJ=-RKNJ 

570  2    RKY=RKY+RKYJ*CP    +RKZJ*SP 

580  RKZ=RKZ-RKYJ*SP    +RKZJ+CP 
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590  RKM=RKM+RKMJ*CP  +RKNJ*SP 

600  RKN=RKN-RKMJ*SP  +RKNJ*CP 

610      1  CONTINUE 

620  YFJ=RKY*TJRC 

630  ZFJ=RKZ*TJRC 

640  YMJ=RKM*TJRC*XJET 

650  ZMJ=RKN*TJRC*XJET 

660  RETURN 

670  END 
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00100 
001  10 
00120+ 


SUBROUTINE  OUTFUT1 

COMMON  YC  9>*DYC  9 > *X*H*WINDXB*UI NDYB*WINDZB* RMACH 

*CPHIV*XFA*YFA*ZFA*THR*XCG*RMASS*PCAN*YCAN*XM*YM*ZM*XXI*YYI 


00130+       *VSND*SPHlW*HTE*GUE*ALPHA*PHlW*XMAX*TBURN*RHO*VEL*rSI*THT* 

00140+  PHI*ZC 14>*DZC 14) 

00150  COMMON    NDIM* I PRI NT*NSTEF*NPRINT* I  COUNT* I  RATE* I  RATI 0*NDIMF 

00160  DIMENSION    OUT( 40* 30 ) *  I OUT< 40* 4) 

00170  100    F0RMAT<F6.3*6F8.2) 

00180  101    F0RMATCF6.3*3F8.2*3F7 .2) 

00190  102  F0RMAT(F6.3*8F7.2) 

00200  103  F0RMATC//*2X**TIME**5X**U**7X**V**7X**W**7X**XI**6X* 

00210+       *YI**6X**HT*> 

00220  104  FORMATC //* 2X* *TIME** 5X* *P** 7X* *Q**7X* *R** 5X* *FSI ** 4X* *THT** 


00230+ 

00240 

00250+ 

00260 

00270 

0G280 

0029J3+ 

00300 

00310 

00320* 

00330 

00340* 

00350* 

00360 

00370 

00380 

00390 

00400 

00410 

00420 

00430 

00440 

00450 

00460 

00470 

00480 

00490  • 

00500 

00510 

00520 

00530 


105 

106 
107 


200 

READ 
10 


4X**FHI*) 

FORMATC //*2X**TIME**2X**ALPHA**3X**PHIW**3X**PHI 1** 

3X**?SI2**2X**  PSI3**2X**  THT4** 4X* *EFB** 3X* *EPBC* ) 

F0RMATCF6.3*4C5X*I 1)) 

FORMATC//** JET  STATES.  1=0N*  0=OFF**/) 
108  FORMAK //* 2X* *TIME** 2X* *WYOKE** 2X* *WCOIL** 4X* *WGY*j 

4X**WGZ**2X**PHI 1D**2X**FSI2D*) 

FORMAT(//**STEP    SIZE=**F6 .4* 2X* *IRATI 0  =  ** I  2) 

DO    10    1  =  1* I  COUNT 

ALL    THE    OUTFUTED   VARIABLES. 

READ<2*>     C0UTCI*J)*J=1*27>*CI0UTCI* J)* J=l*4) 
OUTFUT    THESE    VARIABLES    TO    THE    TTY    IN   APPROPRIATE    COLUMNS    WITH 
HEADINGS. 

PRINT  103 

DO    1     I=1*IC0UNT 

PRINT  100*     (0UT(I*J)*J=1*7) 

PRINT  104 

DO    2    I=1*IC0UNT 

PRINT  101*    OUTCI* 1)*(0UT(I* J)* J=8*13) 

PRINT  105 

DO    3    I=1*IC0UNT 

PRINT  102*    0UT(I*1)*(0UT(I* J)*J=14*21) 

PRINT  108 

DO    5    I=1*IC0UNT 

PRINT  102*    OUTCI* 1)*(0UT(I* J>* J=22*27) 

PRINT  107 

DO    4    I=1*IC0UNT 

PRINT  106*     OUTCI* 1)*(I0UT(I* J)* J=l*4) 

PRINT         200*    H*IRATI0 

STOP 

END 


1 
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00010  SUBROUTINE    INTEG 

00020  COMMON    Y(     9 > * DY (    9 > *X*H*WI NDXB*WI NDYB*WI NDZB*RMACH 

000  30+  *CFHIW*XFA*YFA*ZFA*THR*XCG*RMASS*FCAN*YCAN*XM*YM*ZM*XXI*YYI 

00040+  *VSND*SPHlW*HTE*GUE*ALFHA*PHIW*XMAX*TBURN*RHO*VEL*FSI*THT* 

00050+  PHIj2!04>*DZU4> 

00060  COMMON  NDIM* I FRI NT*NSTEP* NFRINT* I  COUNT* I  RATE* IRATI 0* NDIMF 

00070*  STORAGE  LOCATIONS  FOR  THE  REQUIRED  FAST  STATE  VARIABLES  AND 

00080*  THEIR  DERIVATIVES.  "AUX"  FOR  Y(I)  AND  "BUX"  FOR  Z<I>. 

00090  DIMENSION  AUX  1  (  9>*AUX2(  9>*AUX3(  9>*AUX4(  9)*AUX5(  9)* 

00100+  AUX6C  9>*AUX7(  9)*AUX8(  9 > *BUX 1 ( 1 4) *BUX2( I  4) *BUX3( 1 4) * 

00110+  BUX4( 14>*BUX5( 14>*BUX6( 1 4 > *BUX7 ( 1 4> *BUX8 ( 14) 

00120  DATA  IENTR/l/*IENTRF/2/ 

00140  IF  (IRATE. EG. 1>  GO  TO  100 

00150*  DECIDE  WHETHER  TO  1)  FREDICT  REQUIRED  STARTING  VALUES  BY 

00160*  AN  EULER  BACKSTEP  METHOD*  2)  PERFORM  PREDICTION  OF  THE  Y(I) 

00170*  VARIABLES  DURING  INTEGRATION*  OR  3)  CORRECT  THE  Y(I)  VARIABLES 

00180*  DURING  INTEGRATION. 

00190  IF  (IENTR-2)  20*21*22 

00200*  DETERMINE  REQUIRED  STARTING  VALUES  USING  AN  EULER  BACKSTEF*  AND 

00210*  STORE  INTO  APPROPRIATE  LOCATIONS. 

20220  20  HS=IRATIO*H 

00230  H1=3.*HS 

00240  H3=4./3.*HS 

00250  G1=3.*H 

00260  G3=4./3.*H 

00270  Cl=112./121. 

00280  C2-9./121. 

00290  DO  16  I=1*NDIM 

00300  AUX2(I)=Y(I )-HS*DY(I) 

00310  AUX3( I )=Y(I )-2.*HS*DY(I ) 

00320  AUX4( I )=Y(I )-3.*HS*DY(I ) 

00330  AUX6(I )=DY(I) 

00340  AUX7( I )=DY( I ) 

00350  16  AUX8(I )=0. 

00360  DO  26  I=1*NDIMF 

00370  BUX2CI >=Z(I )-H*DZ(I ) 

00380  BUX3( I )=Z(I >-2.*H*DZ(I> 

00390  BUX4( I )=Z(I )-3.*H*DZ(I > 

00400  BUX6(I)=DZ(I) 

00410  BUX7( I )=DZ(I ) 

00420  26    BUX8( I )=0. 

00430*  PREDICT    THE   Y(I)    VARIABLES. 

00440  21    DO    17    I=1*NDIM 

00450  AUXKI)=Y(I) 

00460  AUX5(I )=DY(I > 

00470  DELT=AUX4( I ) +H3* ( AUX5 ( I ) +AUX5 ( I > -AUX6 ( I >  +AUX7 ( I >  +AUX7 ( I ) > 

00480  Y(I)=DELT-C1*AUX8(I) 

00490  17  AUX8(I)=DELT 

00500  IENTR=3 

00510  GO  TO  100 

00520*  CORRECT  THE  Y(I)  VARIABLES. 

00530  22  DO  18  I=1*NDIM 

00  540  DELT=.125*(9.*AUX1(I)-AUX3(I ) +H 1  * ( DY( I ) +AUX5( I ) +AUX5 ( I ) 

00550+  -AUX6(I>)> 

00560  AUX8(I )=AUX8( I )-DELT 

00570  18    Yd  )=DELT  +  C2*AUX8(I  ) 

00580  DO    19    I=1*NDIM 

00590  AUX7( I )=AUX6(I ) 

00600  •       AUX6( I)=AUX5(I) 
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££610  AUX4CI )=AUX3(I ) 

00620  AUX3CI )=AUX2CI ) 

£0630  19  AUX2CI )=AUX1(I ) 

£064£  IENTR=2 

00650*  NOW  TO  REPEAT  FOR  THE  ZCI)  VARIABLES. 

00660*  DETERMINE  WHETHER  TO  FREDICT  OR  CORRECT  THE  Z(I)  VARIABLES. 

00670  100  IF  (IENTRF.EQ.3)  GO  TO  32 

00680*  FREDICT  THE  ZCI)  VARIABLES. 

00690  DO  27  I=1,NDIMF 

£0700  BUX1CI)=ZCI) 

00710  BUX5CI)=DZCI) 

007  20  DELT=BUX4CI )+G3*(BUX5(I )+BUX5(I )-BUX6CI )+BUX7(I  )+BUX7(I  >) 

00730  ZCI )=DELT-C1*BUX8(I ) 

£0740  27  BUX8CI )=DELT~ 

007  50*  UPDATE    THE   VARIABLE    TIME    BY    "H"    AFTER, EACH    PREDICTION. 
£0760  NSTEF=NSTEF+1 

00770  X=NSTEF*H 

£0780  IPRINT=0 

00790  IENTRF=3 

008  00  RETURN 

008  10*  CORRECT    THE    ZCI)    VARIABLES. 

0£820  32    DO    28    I=1,NDIMF 

008  30  DELT=.125*C9.*BUX1CI )-BUX3CI ) +G 1  * C DZC I ) +BUX5C I )+BUX5C I ) 
00840+  -BUX6CI))) 

££850  BUX8CI )=BUX8CI ) -DELT 

00860  28    ZCI )=DELT+C2*BUX8CI) 

£0870  DO    29    I=1,NDIMF 

00880  BUX7C I )=BUX6CI ) 

£0890  BUX6C I)=BUX5CI ) 

00900  BUX4CI)=BUX3CI ) 

0091E  BUX3CI )=BUX2CI) 

00920  29    BUX2CI )=BUX1CI ) 

009  30  IENTRF=2 
00940  IPRINT=1 
009  50  RETURN 
£0960  END 
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10        SUBROUTINE  I NTER2C C00, C 1 Q, C0 1 , CI \, CB) 

20*  PERFORMS  A  TWO  DIMENSIONAL  INTERPOLATION  BETWEEN  THE  FOUR  CORNERS 

30*  OF  THE  SQUARE  C00-C11  AND  RETURNS  THE  ANSWER  BY  CB . 

40*  COMMUNICATIONS  WITH  SUBROUTINE  AERO  WHICH  CALLS  INTER2. 

50*  DPHI  AND  DALPHA  ARE  THE  INTERPOLATION  INCREMENTS. 

60        COMMON/AEROl /DPHI >  DALPHA 

70        C1=C00+DALPHA*(C10-C00) 

80        C2=C01+DALPHA*(C11-C01) 

90        CB=C1+DFHI*(C2-C1) 

100        RETURN 

110        END 
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READY. 

FORTRAN, OLD, SXDF 

READY. 
RUN,M=1 1000 

03/28/73.     17.05.04 
PROGRAM         SXDF 


INFUT   H,XMAX,NPRINT,IRATIO 

?     .0005*2. ,200, 10 
INPUT   XT0,YT0,HTT0,XTD,YTD,HTTD 

?    9000,3000,100,0,0,0 


TIME 

U 

V 

W 

XI 

YI 

HT 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

.100 

41.41 

0. 

.00 

.00 

0. 

2.07 

.200 

83.14 

1  .29 

1.29 

.02 

.02 

8.29 

.300 

124.61 

-9.71 

-7.74 

.61 

.55 

18.68 

.400 

157.77 

-18.84 

-47.56 

2.57 

1.36 

33.11 

.500 

166.45 

-20.10 

-114.34 

7.15 

2.96 

50.89 

.600 

175.92 

-36.26 

-150.15 

15.35 

5.63 

71.00 

.700 

196.31 

-47.88 

-166.24 

27.98 

9.39 

91.98 

.800 

227.44 

-52.29 

-169.56 

45.23 

14.44 

112.76 

.900 

266.06 

-51  .83 

-162.85 

67.10 

20.97 

132.50 

1.000 

308.77 

-49.27 

-150.62 

9  3.59 

29.04 

150.66 

1.100 

353.54 

-45.26 

-136.15 

124.68 

38.69 

166.91 

1.200 

400.89 

-39.66 

-118.95 

160.36 

49.92 

161  .18 

1.300 

449.40 

-34.79 

-100.35 

200.65 

62.74 

193.53 

1  .400 

498.01 

-30.25 

-86.66 

245.59 

77.18 

204.17 

1.500 

552.44 

-14.23 

-36.64 

295.19 

93.24 

213.61 

1.600 

601.57 

11.12 

37.54 

349.48 

110.87 

222.88 

1.700 

645.79 

27.02 

80.41 

408.31 

129.98 

233.29 

1  .800 

698.37 

16.60 

45.12 

471 .48 

150.43 

245.71 

1.900 

749.20 

-10.08 

-34.29 

539.20 

172.30 

259.58 

2.000 

796.43 

-23.79 

-70.56 

611  .76 

195.74 

273.45 

5h 


TIME 

P 

G 

R 

PSI 

THT 

PHI 

0* 

0* 

0. 

0. 

» 

0. 

» 

90, 

.00 

0. 

i 

•  100 

0. 

0. 

0. 

i 

0. 

t 

90. 

.00 

0. 

•  200 

0. 

-•! 

52 

« 

.52 

45, 

.00 

89. 

.32 

45. 

00 

.300 

0. 

-2.. 

43 

1  . 

.31 

41. 

.34 

77. 

.76 

41  . 

40 

.400 

0. 

-4.; 

24 

""   4 

.16 

18. 

.61 

59. 

.74 

20. 

22 

.500 

0. 

-3.; 

25 

1    4 

.12 

11. 

.62 

35. 

.64 

14. 

61 

.600 

0. 

i 

•l.i 

30 

1. 

.03 

14. 

.15 

20. 

.65 

1.5. 

77 

.700 

0. 

-l.i 

30 

4 

.64 

16. 

17 

11. 

.83 

16. 

32 

.800 

0. 

-•! 

53 

4 

.19 

17. 

.03 

6. 

.90 

16. 

46 

.900 

0. 

T  .  '< 

25 

4 

.23 

17. 

59 

44 

.18 

16. 

52 

1  .000 

0. 

-.23 

4 

.03 

17. 

.90 

2. 

.58 

16. 

54 

1.100 

-     0. 

•  . 

15 

4 

.11 

18. 

13 

1. 

.45 

16. 

55 

1.200 

0. 

. 

11 

< 

.09 

18. 

.23 

1. 

.12 

16. 

55 

1.300 

0. 

—  . 

10 

4 

.07 

18. 

.40 

1. 

.03 

16. 

.55 

1.400 

0. 

. ! 

26 

"*   4 

.06 

18. 

.39 

4 

.92 

16. 

.55 

1.500 

0. 

1  . 

16 

•"   4 

.39 

18. 

.30 

5. 

.77 

16. 

.54 

1.600 

0. 

1  . 

15 

"    < 

.40 

17. 

.99 

13. 

.22 

16, 

.49 

1.700 

0. 

. 

17 

—    4 

.09 

174 

.68 

17, 

.58 

16. 

.41 

1.800 

0. 

-.« 

38 

•  31 

17. 

.66 

14, 

.66 

16. 

.40 

1.900 

0. 

-1  .1 

37 

4 

.37 

17. 

.90 

7< 

.83 

16. 

.45 

2.000 

0. 

•  i 

32 

4 

.01 

18, 

.12 

4. 

.35 

16. 

.47 

TIME 

ALPHA 

FHIW 

PHI1 

PSI2 

PSI3 

THT4 

EPB 

EPBC 

0. 

0* 

0 

• 

71 

.57 

89 

.40 

0 

. 

0 

t. 

•  00 

•  00 

•  100 

0. 

0 

. 

71 

.54 

89 

.28 

.13 

•  00 

.03 

•  13 

•  200 

1.25 

45 

.00 

71 

.57 

88 

•  88 

- 

.02 

.04 

.05 

•  04 

.300 

5.69- 

•128 

.55 

71 

.86 

78 

.33 

- 

.10 

.17 

.03 

.18 

.400 

17.96- 

•158 

.39 

70 

.92 

59 

.40 

- 

•  10 

•  80 

.03 

.81 

.500 

34.90- 

•170 

.03 

64 

.24 

36 

.08 

- 

•  16 

.32 

•  02 

•  35 

.600 

41.28- 

•166 

.42 

62 

.46 

21 

.07 

- 

.18 

.09 

•  02 

.21 

.700 

41.39- 

■163 

.93 

62 

.35 

12 

.18 

- 

.19 

.10 

.02 

.20 

.800 

37.96- 

■162 

.86 

62 

.43 

7 

.30 

- 

.20 

.02 

•  03 

.19 

.900 

32.71- 

■162 

.35 

61 

.75 

4 

.64 

- 

.16 

•  04 

.03 

.19 

1.000 

27.17- 

■161 

.89 

62 

.99 

3 

.10 

- 

.17 

•  02 

•  01 

.16 

1.100 

22.09- 

■161 

.61 

62 

.93 

2 

.03 

- 

•  16 

.02 

.01 

.15 

1.200 

17.37- 

•161 

.56 

62 

.97 

1 

.59 

.06 

.05 

.04 

.08 

1.300 

13.30- 

■160 

.88 

63 

.01 

1 

.67 

- 

.10 

•.26 

.03 

.25 

1.400 

10.44- 

•160 

.75 

63 

.05 

1 

.51 

•  06 

.19 

.03 

.22 

1.500 

4»07- 

■158 

.77 

63 

.43 

6 

•  23 

•  20 

.97 

.02 

.98 

1.600 

3.72 

16 

.50 

71 

.49 

13 

.76 

.23 

.02 

.03 

•  24 

1»700 

7.48 

18 

.57 

71 

.23 

18 

.22 

.21 

.02 

•  02 

.22 

1.800 

3.94 

20 

•  20 

70 

.80 

15 

.76 

.16 

•  09 

•  02 

.17 

1.900 

2.73- 

■163 

.62 

69 

.92 

9 

•  06 

- 

.17 

•  04 

.03 

.19 

2.000 

5.34- 

•161 

.37 

69 

.50 

5 

.64 

.15 

•  00 

.03 

•  15 
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TIME 

WYOKE 

WCOIL 

WGY 

WGZ 

PHI  ID 

PSI2D 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

•  100 

-.02 

.03 

.02 

•  03 

-.02 

•  03 

.200 

.00 

.19 

-.06 

-.01 

•  00 

-.48 

.300 

.10 

-.02 

.01 

-.02 

.10 

-2.73 

.400 

-.81 

.06 

-.03 

.02 

-.81 

-3.90 

.500 

-.70 

.02 

.01 

.04 

-.70 

-3.39 

.600 

-.03 

.03 

.01 

.03 

-.03 

-2.04 

.700 

•  20 

-.01 

•  01 

-.01 

.20 

-1.20 

.800 

-.03 

-.01 

-.04 

-.02 

-.03 

-.57 

.900 

•  21 

.06 

-.10 

.04 

•  21 

-.27 

1.000 

-.16 

.03 

-.09 

.00 

-.16 

-.19 

1.100 

-.01 

•  00 

.01 

.01 

-.01 

-.18 

1  .200 

-.01 

-.04 

-.04 

-.01 

-.01 

.02 

1.300 

-.01 

-.01 

.02 

.01 

-.01 

-.13 

1.400 

-.01 

-.07 

-.04 

.03 

-.01 

.20 

1.500 

1.60 

.00 

-.01 

-.03 

1  .60 

1.22 

1.600 

-.08 

-.02 

-.03 

.05 

-.08 

1  .20 

1.700 

-.03 

-.00 

-.04 

-.02 

-.03 

•  19 

1.800 

-.07 

.02 

-.06 

•  01 

-.07 

-1  .01 

1.900 

-.54 

.05 

-.03 

.05 

-.54 

-1.08 

2.000 

-.14 

-.07 

-.03 

.05 

-.14 

-.05 

JET    STATES.     1=0N,    0=OFF 


0. 

0 

0 

0 

0 

.100 

0 

0 

0 

0 

.200 

0 

0 

1 

.300 

0 

1 

0 

•  400 

0 

0 

1 

.500 

1 

0 

0 

1 

.600 

0 

1 

0 

.700 

0 

1 

0 

•  800 

0 

0 

1 

.900 

0 

1 

0 

1.000 

0 

0 

1 

1.100 

0 

1 

0 

1.200 

0 

1 

0 

1.300 

0 

1 

0 

1.400 

0 

0 

0 

0 

1.500 

0 

0 

0 

0 

1.600 

0 

0 

0 

0 

1.700 

0 

0 

0 

0 

1.800 

0 

0 

0 

0 

1.900 

0 

0 

0 

0 

2.000 

0 

0 

0 

0 

STEP  SIZE= 

.0005 

IRATIO=10 

STOP. 
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APPENDIX  B 
PROGRAM  VARIABLES 


The  following  is  an  alphabetical  listing  of  the  computer  program  variables 

except  for  the  variables  in  subroutine  SEEKER. 

Missile  angle  of  attack 

Angle  of  attack  entry  argument  for  coefficient  table  lookup 

Storage  locations  for  integration  routine  used  for  HS. 


ALPHA 

ALPHAE 

AUX  1 

• 

AUX  8 

BUX  1 

• 

BUX  8 

CI,  C2 

COO,  CIO 

C01,  Cll 

CB 

CL 

CM 

CMP 

CMPB 

CN 

CNP 

CNPB 

C0N1 

conU 

CP 

CPE 

CPHI1 

CPHTW 

CX 

CXB 

CY 


Y.,  Y.  - 
i'  l-l 
respectively 


■  Yi-2'  Yi-3' 


DY.,  DY.  _, 
l    l-l' 


dy 


i-2 


truncation  corrector 


Storage  locations  for  integration  routine  used  for  H. 

V  zi-i'  h-2>  h-y  DZ 

respectively 


,  DZ.  _,  DZ.  _,  truncation  corrector 
l-l    i-2 


Intermediate  answers  for  2-D  coefficient  lookup 
Entry  points  for  2-D  coefficient  lookup 

Coefficient  value  return  by  2-D  coefficient  lookup 

Aerodynamic  total  rolling  moment  coefficient  in  body  fixed  axes 

Aerodynamic  total  pitching  moment  coefficient  in  body  fixed  axes 

Aerodynamic  total  pitching  moment  coefficient  in  pitch  axes 

Array  of  pitching  moment  coefficient  for  baseline  missile  (no 
JRC)  vs  a  and  $ 

Aerodynamic  total  yawing  moment  coefficient  in  body  fixed  axes 

Aerodynamic  total  yawing  moment  coefficient  in  pitch  axes 

Array  of  yawing  moment  coefficient  for  baseline  missile  (no 

JRC)  vs  a  and  $ 
w 

Program  constants 


Dummy  variable 

Vector  of  cosines  of  1^  for  each  engine 

Cos  (O 

Cos  (fy) 

Aerodynamic  total  X  force  coefficient  in  body  fixed  axes 

Array  of  X  force  coefficient  for  baseline  missile  (no  JRC) 

vs  a  and  $ 
w 

Aerodynamic  total  Y  force  coefficient  in  pitch  axes 
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CYP 
CYPB 

cz 

CZP 

CZPB 

DALPHA 

DELT 

DIA 

DPHI 

DY 

DZ 

EO 


Aerodynamic  total  Y  force  coefficient  in  "body  fixed  axes 

Array  of  Y  force  coefficients  for  baseline  missile  (no  JRC) 

vs  a  and  $ 
w 

Aerodynamic  total  Z  force  coefficient  in  body  fixed  axes 
Aerodynamic  total  Z  force  coefficient  in  pitch 

Array  of  Z  force  coefficients  for  baseline  missile  (no  JRC) 

vs  a  and  $ 
w 

The  residule  of  a    required  for  2-D  coefficient  lookup 

Intermediate  variable  in  integration  routine 

Missile  diameter 

The  residule  of  $   required  for  2-D  coefficient  lookup 

Vector  of  derivatives  of  Y 

Vector  of  derivatives  of  Z 

Quaterion  variables 


E3 

EOD 


Derivatives  of  E0-E3 


E3D 

EPB 

EPBC 

EPBCO 

EPS 

Gl,  G3 

GRAV 

HI,  H3 

H 

HS 

HT 

HTD 

HTMT 

HTTO 

HTTD 

HTTI 

IA 


Angular  error  between  gyro  spin  axis  and  target  line  of  sight 
in  radians 

Angular  error  between  coil  housing  axis  and  target  line  of 
sight  in  radians 

EPBC  in  degrees  for  outputting 

Error  due  to  non-orthogonality  of  E0-E3 

Internal  constants  for  integration  routine 

Gravity  in  ft /sec 

Internal  constants  for  integration  routine 

Step  size  for  fast  integration  routine 

Step  size  for  slow  integration  routine 

Altitude  of  missile 

Derivative  of  HT 

Altitude  difference  between  missile  and  target 

Initial  height  of  target 

Target  climb  rate 

Altitude  of  target 

Entry  location  on  a  for  2-D  coefficient  lookup 
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ICOUNT 
IENTR 

IENTRF 

IJET 

IOUT 

IP 

IPRINT 

IRATE 

IRATIO 

IRT2 

NDIM 

NDIMF 

NPRINT 

NSTEP 

OUT 

P 

PCAN 

PD 

PHI 

PHI1 

PHTIO 

PHI1D 

PHIE1 

PHIJ 

PHIW 

PHIWT 

PHIWTE 

PRMT 

PRMTY 

PRMTZ 

PSI 

PS  12 

PS  120 

PSI2D 

PS  13 


Number  of  times  outputting  is  performed 

1 — Calculate  required  starting  values,  2 — Predict  Y(l)  for 
integration,  3 — Correct  Y(l)  for  integration 

2 — Predict  Z(l)  for  integration,  3 — Correct  z(l)  for  integration 

Array  for  each  jet;  0--off,  1 — on 

Array  for  outputting  interger  variables 

Entry  location  on  $  for  2-D  coefficient  lookup 

0 — No  outputting  permitted,  1 — Outputting  permitted 

Controls  integration  routine;  0 — Integrate  both  Y  and  Z,  1 — Integrate 
only  Z 

HS/H 

2*   IRATIO 

Dimension  of  Y(l) 

Dimension  of  Z(l) 

Outputting  interval  =  NPRINT*H 

Counter  for  number  of  H  integration  steps 

Array  for  outputting  real  variables 

Missile  roll  rate  in  body  axes  in  rad/sec 

Pitch  canard  deflection 

Derivative  of  P 

$  in  degrees 

§_  in  radians 

PHI1  in  degrees  for  outputting 

Derivative  of  PHI1 

3L  for  thrust er  number  1 

$  for  thruster  in  degrees 
J 

$  in  degrees 

$  reduced  to  a  range  of  0  -  9°  degrees 

$  entry  argument  for  coefficient  table  lookup 

Dummy  variable 

Control  function  for  JRC's  aligned  with  missile  Y  axis 

Control  function  for  JRC's  aligned  with  missile  Z  axis 

Y  in  degrees 

Y?  in  radians 

Yp  in  degrees  for  outputting 

Derivative  of  PS 12 

y     in  radians 
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PS  130 
P3I3D 

Q 

QD 

QUE 

R 

RD 

RHO 

RKM 

RKMB 

RKN 

RKNB 

RKY 

RKYB 

RKZ 

RKZB 

RMACH 

RMASS 

RMASSO 

RMASSD 

S 

SIGA 

SIGAD 

SIGB 

SIGBD 

SP 

SFE 

SPHI1 

SPHLW 

Til 


Y„  in  degrees  for  outputting 

Derivative  of  P3I3 

Missile  pitch  rate  in  body  axes  in  rad/sec 

Derivative  of  Q, 

Dynamic  pressure  ^pv  in  lbs /ft 

Missile  yaw  rate  in  body  axes  in  rad/sec 

Derivative  of  R 

p  in  slugs /ft 

Total  moment  amplification  factor  along  missile  Y  axis 

Array  of  moment  amplification  factor  vs  a  and  $ 

J 

Total  moment  amplification  factor  along  missile  Z  axis 
Array  of  moment  amplification  factor  vs  a  and  $T 
Total  force  amplification  factor  along  missile  Y  axis 
Array  of  forces  amplification  factor  vs  a  and  $T 
Total  force  amplification  factor  along  missile  Z  axis 
Array  of  forces  amplification  factor  vs  a  and  $T 
Mach  number  of  missile 
Instantaneous  missile  mass 
Initial  missile  mass 
Derivative  or  RMASS 


rrd 

"IT 


Missile  reference  area  = 

a  in  degrees 

Derivative  of  SIGA 

ov.  in  degrees 

Derivative  of  SIGB 

Dummy  variable 

Vector  of  sines  of  $^  for  each  engine 

Sin  ($1) 

Sin  ($w) 

Elements  of  matrix  transformation  from  body  to  inertial  axes 


T33 

TBJRC 

TBURN 

THT 

THTU 


Burn  time  of  JRC's 

Burn  time  of  missile  main  thruster 

9  in  degrees 

Q*    in  radians 
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THR 

THT^O 

THT^D 

TJRC 

TSJRC 

U 

UD 

UT 

V 

VD 

VEL 

VEL1 

VEL2 

VSND 

VT 

W 

wcz 

WGY 
WGZ 
WINIKB 

WINDZB 
WINDXI 

• 

WUTDZI 

WT 

X 

XCG 

XCGO 

XCGD 

XF 

XFA 

XFB 

XI 

XID 

XIMT 

XJET 


Main  thrust er  thrust 

0^  in  degrees  for  outputting 

Derivative  of  THTU 

Thrust  of  one  JRC 

Enable  time  for  JRC ' s 

U  velocity  of  missile 

Derivative  of  U 

UT 

V  velocity  of  missile 

Derivative  of  V 
Total  missile  velocity 
Dummy  variable 
Dummy  variable 
Velocity  of  sound 

VT 

W  velocity  of  missile 

Q  of  seeker  cage  along  Z  axis 

Q  of  seeker  gyro  along  Y  axis 

Q  of  seeker  gyro  along  Z  axis 

Components  of  wind  along  missile  body  axes 

Components  of  wind  along  inertial  axes 


WT 

Time 

Position  of  missile  C.G.  relative  to  aerodynamic  reference  point 

Initial  XCG 

Derivative  of  XCG 

Total  force  along  X  axis 

Total  aerodynamic  force  along  X  axis 

Total  aerodynamic  force  along  X  axis  for  baseline  missile  (no  JRC 

Inertial  X  position  of  missile 

Derivative  of  XI 

Inertial  X  position  of  target  relative  to  missile 

Position  of  JRC's  from  aerodynamic  reference  point 
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XM  Total  moment  along  X  axis 

XMAX  Total  run  time 

XMB  Total  moment  along  X  axis  for  baseline  missile  (no  JRC) 

XMT  X  position  of  target  relative  to  missile  in  body  axes 

XMTD  Derivative  of  XMT 

XTO  Initial  XTI 

XTD  Derivative  of  XTI 

XTI  Inertial  X  position  of  target 


XXI 

XXIO  Initial  3 


XXID  Derivative  of  XXI 

Y  Vector  of  slow  integration  state  variables 

YCAN  Deflection  of  yaw  canard 

YF  Total  along  Y  axis 

YFA  Total  aerodynamic  force  along  Y  axis 

YFB  Total  aerodynamic  force  along  Y  axis  for  baseline  missile  (no  JRC) 

YFJ  Total  force  along  Y  axis  due  to  JRC's 

YI  Inertial  Y  position  of  missile 

YID  Derivative  of  YI 

YIMT  Inertial  Y  position  of  target  relative  to  missile 

YM  Total  moment  along  Y  axis 

YMB  Total  moment  along  Y  axis  for  baseline  missile  (no  JRC) 

YMJ  Total  moment  along  Y  axis  due  to  JRC's 

YMT  Y  position  of  target  relative  to  missile  in  body  axes 

YMTD  Derivative  of  YMT 

YTO  Initial  YTI 

YTD  Derivative  of  YTI 

YTI  Inertial  Y  position  of  target 

YYI  Iyy 

YYIO  Initial  YYI 

YYID  Derivative  of  YYI 

Z  Vector  of  fast  integration  state  variables 

ZF  Total  along  Z  axis 

ZFA  Total  aerodynamic  force  along  Z  axis 
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ZFB  Total  aerodynamic  force  along  Z  axis  for  baseline  missile  (no  JRC) 

ZFJ  Total  force  along  Z  axis  due  to  JRC ' s 

ZM  Total  moment  along  Z  axis 

ZMB  Total  moment  along  Z  axis  for  baseline  missile  (no  JRC) 

ZMJ  Total  moment  along  Z  axis  due  to  JRC's 

ZMT  Z  position  of  target  relative  to  missile  in  body  axes 

ZMTD  Derivative  of  ZMT 
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